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ABSTRACT 


A unified  approach  for  the  implementation  of  direct  time  integration 
procedures  in  structural  dynamics  is  presented.  Two  key  performance 
assessment  factors  are  considered,  viz.  , computational  effort  and 
error  propagation.  It  is  shown  that  these  factors  are  strongly  affected 
by  details  in  the  reduction  of  the  second-order  equations  of  motion  to  a 
system  of  first-order  equations,  and  by  the  computational  path  followed 
at  each  time  step.  Part  I is  primarily  devoted  to  the  study  of  the 
organization  of  the  computational  process.  Specific  implementation 
forms  derivable  from  the  unified  approach  are  studied  in  detail,  and 
rated  accordingly. 
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Section  1 
INTRODUCTION 


Motivation 

The  use  of  direct  time  integration  procedures  for  the  numerical  treat- 
ment of  structural  dynamics  problems  has  attracted  an  increasing  number 
of  investigators  in  recent  years.  Most  studies  have  concentrated  on  the 
development  of  new  methods,  and  on  the  formulation  and  verification  of 
performance-evaluation  criteria.  These  studies  have  established  sta- 
bility and  accuracy  (in  that  order)  as  the  most  important  attributes  of  a 
time  integration  operator,  and  there  is  an  emerging  consensus  on  the 
use  of  evaluation  criteria  appropriate  to  the  treatment  of  mechanical 
equations  of  motion  (or,  in  fact,  any  stiff-oscillatory  system).  Accuracy 
can  be  assessed  by  measuring  numerical  damping  and  frequency  distor- 
tion as  a function  of  the  sampling  rate  when  the  integrator  is  applied  to 
undamped  linear -oscillator  models.  Stability  can  be  assessed  by  studying 
the  amplitude  growth  of  the  computed  solutions  of  linear  and  nonlinear 
model  problems. 

There  is  presently  widespread  agreement  on  the  stability  advantages  of 
implicit  operators  for  treating  structural  dynamics  problems,  in  which 


the  gross  response  of  the  discrete  model  is  dominated  by  the  low-frequency 
components  of  the  modal  spectrum.  Explicit  schemes  are  still  widely  used, 
however,  in  problems  involving  shock-type  excitation  and  wave-propagation 
type  response,  where  the  contribution  of  the  high-frequency  components 
remain  significant  throughout  the  analysis  process.  Performance  cross- 
over points  between  explicit  and  implicit  schemes  have  not  been  firmly 
established  as  yet. 


The  algorithmic  properties  pertaining  to  stability  and  accuracy  play  an 
important  role  in  the  preliminary  evaluation  of  time  integration  operators, 
and  represent  a useful  source  of  information  for  selecting  methods  appro- 
priate to  different  classes  of  problems.  There  have  been  conflicting 


reports,  however,  as  regards  the  actual  performance  of  the  same 
integrator  in  different  computer  programs.  Contradictory  claims  are 
especially  frequent  in  nonlinear  analysis*  These  unresolved  questions 
suggest  that  the  computer  implementation  of  a time  integration  procedure 
has  a substantial  bearing  on  the  performance  of  a transient  analysis  pro- 
gram, It  will  be  shown  that  computer  implementation  details  affect  two 
important  performance  characteristics:  the  computational  efficiency  and 
the  propagation  of  computational  error. 


Computational  Efficienc* 


For  a given  stepsize  sequence,  the  efficiency  of  a transient  analyzer  is 
largely  controlled  by  the  volume  of  arithmetic  operations  and  data  transfers 
performed  at  each  step.  These  activities  are  governed  in  turn  by  the  allo- 
cation of  storage  resources  and  by  the  order  and  manner  in  which  array- 
processing operations  are  carried  out.  It  will  be  shown  that  several  ways 
of  implementing  the  basic  time -advancing  cycle  are  possible  even  if 
the  same  integration  operator  is  used,  and  that  these  implementations 
display  varying  efficiency  characteristics. 

Computational  Error  Propagation 


Computational  errors  arise  from  three  main  sources: 

• Rounding  errors  resulting  from  the  use  of  finite -precision  arith- 
metic. This  source  is  significant  if  the  linear  or  nonlinear 
algebraic  equations  that  must  be  solved  at  each  time  station  are 
ill-conditioned. 

• Inaccuracies  resulting  from  the  approximate  solution  of  the 
aforementioned  equations  by  iterative  techniques  in  the  case  of 
nonlinear  problems. 


Errors  in  the  calculation  of  applied  forces.  This  source  can  be 
significant  if  the  excitation  depends  on  the  solution  of  a coupled 
initial  value  problem,  as  in  the  case  of  fluid-structure  interaction, 
dynamic  thermal  loads,  and  certain  types  of  material  nonlinearity. 
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The  "local"  computational  error  due  to  the  combined  effect  of  the  pre- 
ceding sources  at  a typical  time  station  t^  may  be  viewed  as  a spurious 
excitation  source.  This  disturbance  is  propagated,  and  usually  amplified, 
by  the  integration  operator  to  the  stream  of  computed  solutions  for  t 
The  effect  of  a local  perturbation  on  downstream  solutions  will  be  called 
the  propagated  computational  error.  The  overall  computational  error  at 
a given  time  station  results  from  the  combination  of  the  local  error  and 
of  the  propagated  errors  emanating  from  all  previous  stations;  this 
quantity  will  be  called  the  accumulated  (or  global)  computational  error. 

A comprehensive  theory  of  error  propagation  in  the  direct  integration  of 
nonstiff  ordinary  differential  equations  (ODE)  has  been  laid  out  by  Henrici 
[1,  2]  . His  investigations  concentrated  on  the  development  of  asymptotic 
bounds  for  the  accumulated  computational  error,  i.  e.  , bounds  applicable 
if  lu).  hl«  1,  where  h is  the  integration  stepsize  and  are  characteristic 
roots  of  the  associated  homogeneous  system.  This  paper  complements 
and  extends  Henrici' s results  to  the  stiff  ODE  systems  arising  in  struc- 
tural dynamics  applications,  and  studies  the  crucial  effect  of  implementa- 
tion details  on  error  propagation. 

Objectives  of  the  Paper 


I i The  paper  has  two  basic  objectives.  First,  to  present  a general  framework 

for  reducing  the  second-order  equations  of  structural  dynamics  to  an  equi- 
! valent  first-order  system,  to  which  a pair  of  linear  multistep  (LMS) 

( implicit  time  integrators  is  applied,  and  to  categorize  the  resulting 

t I 

I j implementation  variants.  This  approach  is  shown  to  produce  many  of  the 

i computer  implementations  in  current  favor,  as  well  as  a myriad  of 

' hitherto  untried  ones.  Second,  to  study  the  dependence  of  the  propagated 

computational  error  on  various  implementations  resulting  from  the  general 
approach,  and  to  correlate  that  behavior  with  compositional  properties  of 
the  integration  operator. 


Because  of  its  length  and  varying  degree  of  mathematical  sophistication 
expected  from  the  reader,  the  paper  is  subdivided  into  two  parts,  which 
address  the  preceding  objectives  in  turn. 
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Outline  of  Part  I 

The  reduction  of  the  governing  matrix  equations  of  linear  structural 
dynamics  to  a first-order  system  is  presented  following  the  auxiliary- 
vector  procedure  introduced  by  Jensen  [3].  Two  LMS  implicit  integration 
operators  are  applied  to  the  resulting  first-order  system  to  derive  a set 
of  coupled  difference  equations.  This  set  contains  arbitrary  weighting 
matrices  associated  with  the  choice  of  auxiliary  vector.  The  selection  of 
convenient  weighting  matrices  leads  to  specific  formulations,  which  in  turn 
subdivide  into  variants  according  to  the  computational  path  followed  in  the 
typical  time -advancing  cycle.  Operation  counts  associated  with  major  array - 
manipulation  tasks  are  exhibited  as  a function  of  the  computer  implementa- 
tion, sparseness  attributes  of  the  mass  and  damping  matrices,  and  certain 
compositional  properties  of  the  integration  operators.  Associated  starting 
procedures  that  minimize  the  need  for  additional  computations  are  des- 
cribed. Finally,  the  advantages  and  disadvantages  of  the  various  formulations 
are  summarized. 

Practical  considerations  on  the  computational  error  propagation  behavior 
associated  with  various  implementations  are  interspersed  throughout  Part  I. 
These  considerations  represent  an  extract  of  relevant  conclusions  of  the 
investigations  reported  in  Part  II,  and  are  offered  so  as  to  make  Part  I 
relatively  self-contained  for  those  readers  that  are  not  particularly  moti- 
vated to  follow  the  details  of  the  error  analysis. 

Linear  dynamic  systems  are  emphasized  throughout,  as  this  keeps  the 
major  points  of  the  exposition  uncluttered  from  unnecessary  complexities 
brought  about  by  the  intrusion  of  nonlinear  terms.  Many  of  the  conclusions 
regarding  computer  implementation  aspects  and  error  propagation  behavior 
are  applicable  ''locally"  to  the  linearized  forms  of  the  governing  nonlinear 
dynamic  equations,  whether  resulting  from  the  tangent-stiffness  or  the 
pseudo-force  approach. 
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Section  2 

SUMMARY  OF  NEW  RESULTS  IN  PART  I 


Jr  r 


1.  Categorization  of  computer  implementation  procedures  according 
to  the  selection  of  auxiliary  vector  for  the  reduction  to  a first- 
order  system. 

2.  Distinction  of  three  basic  computational  paths  for  the  evaluation 
of  auxiliary  vectors  in  the  time -advancing  cycle. 

3.  Tabulation  of  computational  effort  associated  with  various 
implem  entation  s . 

4.  The  concept  of  '^starting  operator  family'*  to  minimize  the  need 
for  additional  computations  in  the  starting  procedure. 


5.  Formal  embedding  of  the  special,  two-derivative  operators  of 

structural  mechanics  (Newmark,  Houbolt,  Wilson)  within  the 
general  implementation-classification  approach. 
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Section  3 

FORMULATION  OF  DIFFERENCE  EQUATIONS 


Governing  Equations 

The  governing  equations  of  motion  of  a linear,  discrete  mechanical 
system  with  n^  degrees  of  freedom  may  be  written 

Mu  + Du+Ku  = £{t) 

where  u(t)  is  a state  vector  containing  n^  displacement  coordinates 
defining  the  spatial  configuration  of  the  discrete  system  at  each  time  t 
M,  D,  and  K are  the  linear  mass,  damping,  and  stiffness  matrices, 
respectively,  and  £(t)  is  a vector  of  time-varying  applied  generalized 
forces:  in  (1)  a dot  denotes  temporal  differentiation. 

Reduction  to  a First-Order  System 

A general  procedure  for  reducing  (1)  to  a system  of  first-order  equations 
has  been  presented  by  Jensen  [3].  He  introduced  an  auxiliary  vector  v 
defined  as  a weighted  combination  of  displacements  and  velocities: 

v = AMu  + Bu  (2 


where  A and  B are  (for  the  moment)  arbitrary  n.  by  n^  matrices  with  the 
^ ^ I t 

proviso  that  A be  nonsingular.  The  resulting  first-order  system  is 


AM  0 M u 


A D - B I 


A K 0 


where  I and  0 denote  the  identity  and  null  matrix,  respectively.  Eq.  (3) 
fits  within  the  standard  ODE -system  format: 


i = F(^,t) 


(4) 


where 


A 1 
A D 

r*-' 


- B I 


V - B u 


A (f  - K u) 


Integration  Method 

We  consider  the  numerical  integration  of  the  first-order  system  (3)  using 
m-step,  one -derivative,  linear  multistep  (LMS)  operators  to  discretize 
the  time  variation  of  u and  v.  A detailed  treatment  of  these  methods  can 
be  found  elsewhere.  I g. . [1-2,  4-6].  For  a constant  stepsize  h.  LMS 
operators  applied  to  Eq.  (3)  at  the  n-th  time  station  t^  can  be  written 


—A 


V.  V 

H -n-i 


Hi 


^ 1,  n S m 


O'  0 

O o 


where  cy.  . P-  , Y • ^nd  6^  are  scalars  associated  with  specific  finite 

^ . j , ,T  stand  for  the  vectors  n(t, 

difference  operators,  andu^,  stana  lu 

computed  through  the  discretization  procedure  (7).  The  above 

mffe’rence  expressions  may  be  scaled  by  arbitrary  factors.  In  the  sequel 

they  will  be  normalized  by  requiring  that 

_ , (8) 

= '^o  - ^ 

Remark.  In  the  most  general  case  one  could  use  an  m-step  and  a k-step 
(k  < m)  LMS  operator  in  (7).  The  latter  can  always  be  considered  as  an 
m-3tep  operator,  however,  by  appropriately  appending  (m-k)  zero  coeffi- 

cients. 
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With  a view  to  subsequent  manipulations,  it  is  convenient  to  introduce  a 

more  compact  notation  than  (7),  which  nevertheless  separates  explicitly 

the  current  state  variables  u , u , v , v from  historical  terms 
■ “n  ^n  “n  ^n  ■ 

involving  past  solutions: 


u 

— n 

+ 

— n 

= h u + 
P n 

h b"" 

— n 

V 

— n 

+ 

— n 

= h,  V + 

6— n 

hd^ 
— n 

generalized 

stepsizes" 

hp.h 

0 

CD. 

U 

^6  = 

h 6 

o 

and  vectors  a,  b,  ...  denote  linear  combinations  of  past  solutions: 


z , z z ,z 
abed 
— n --n  — n * -n 


^'i  Pi  *1 


o 6 Y 6 
m m 


in  which  z is  an  arbitrary  vector  symbol,  which  may  stand  for  u , . 

etc.  The  current  state  solution  (u  , v ) can  be  expressed  from  (9)  in 

— n — n 

the  form 


h^  a 


where 


will  be  called,  following  Jensen  [3]  , the  historical  vectors  pertaining  to 
u and  V,  respectively. 
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Implicit  Discretization 


We  shall  consider  only  implicit  LMS  operators,  for  which  both  and  6^ 

are  nonzero,  whence  h^  ^ 0 and  h.  ^ 0.  Elimination  of  a and  v from  (3) 

p 0 “H.  — n 

and  (12)  then  produces  the  linear  algebraic  system  of  order  2n^: 


A M + h^  B -h^  ; 

rw  ^ P ^ P r 

AD-B+hpAK  TII 


A M h“  I 

— n I 

(A  D - B)  + l\h^  + h„  A f ( 
^ — n n p ^ “Ti  I 


where  'n=h„/h=P  /6  . Finally,  elimination  of  v from  (14)  yields  the 
P 0 o o —n 

n^-th  order  linear  system 


where 


E = M + D + h_  h K + (h  - h ) C 
e = [M  + h JD-C)]  + h-A'^  h'"  + h- h f 


A’^  B 


As  M,  D,  and  K are  normally  symmetric,  C = A ^ B should  also  be  a 

/-s_/  ^s/ 


/“S-/  r 

symmetric  matrix  if  h i h in  order  to  avoid  the  appearance  of  an 

P 0 

unsymmetric  coefficient  matrix  E.  This  may  pose  certain  constraints  on 
the  selection  of  A and  B. 


The  generality  of  the  discretization  (7),  in  which  different  LMS  operators 
may  be  used  for  u and  v,  has  two  important  advantages.  First,  the  use  of 
different  approximations  for  u and  v is  allowed,  should  that  prove  desirable 
on  account  of  the  physical  meaning  of  v resulting  from  specific  choices  of 
A and  B in  (2).  Second,  special  two-derivative  integrators  devised  for 
treating  the  second-order  system  (1)  directly,  such  as  the  widely  used 
Newmark  [6],  Houbolt  [7],  and  Wilson  [8 ] operators , can  be  presented 
as  special  cases  of  a slight  extension  of  (7),  in  which  a historical  u-term 
is  appended  to  (7b),  as  shown  in  Appendix  A.  It  follows  that  a separate  study 
of  the  computer  implementation  of  those  special  methods  is  not  required. 
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The  selection  of  a common  LMS  operator  for  both  u and  v is  natural  if 
the  analyst  wishes  to  keep  similar  discretization  accuracy  on  both  state 
vectors.  For  this  particular  choice,  Eqs.  (9)  become 
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Section  4 

COMPUTEE  IMPLEMENTATION 

Implementations  of  the  implicit  time  integration  procedure  outlined  in 
the  preceding  section  may  differ  in  two  respects:  (a)  the  selection  of 

weighting  matrices  A and  B in  (2),  and  (b)  the  manner  in  which  the  result- 
ing  auxiliary  vector  v and  its  temporal  derivative  v are  computed  at  each 
step.  We  deal  with  the  latter  topic  first  before  passing  to  discuss  the 
selection  of  computationally  convenient  A and  B. 

Computational  Paths 


Three  basic  computational  paths,  labelled  (0),  (1)  and  (2),  may  be  followed 
in  advancing  the  discrete  solution  over  a typical  time  step.  The  associated 
sequences  of  calculations,  for  arbitrary  A and  B,  are  flow  charted  in 
Figure  1. 

The  path-identification  index  (0-2)  gives  an  idea  of  the  number  of  backward 
difference  operations  performed  in  the  determination  of  andjy^  over  each 
solution-advancing  cycle.  As  shown  in  Part  II,  this  index  plays  an  impor- 
tant role  in  the  study  of  the  error  propagation  behavior  of  particular 
implementations. 

Path  (0)  is  consistent  with  the  first-order  difference  system  (14),  i.  e. , 
the  computed  vectors  satisfy  (14)  exactly  if  computational  errors  are 
ignored.  The  original  ODE  system  (1)  is  also  satisfied.  On  the  other  hand, 
the  computed  u^  , ^ and  do  not,  in  general,  verify  the  auxiliary -vector 
definition  (2),  which  is  satisfied  only  in  the  limit  h-*  0, 

A slight  variant  of  path  (0),  labelled  (0*),is  also  shown  in  Figure  1.  In 
this  variant  the  velocity  vector  ^ ^ is  recomputed  so  that  (2)  is  also 
verified  at  past  stations.  This  additional  operation  is  not  only  expensive 
for  general  A and  B,  but  worthless  as  it  does  not  improve  the  error  propa- 
gation characteristics  of  the  integrator.  It  is  included  here,  however, 
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STARTING  PROCEDURE: 


Computational  paths  for  arbitrary  auxiliary  vector 


because  it  occurs  naturally  in  the  conventional  choice  v = u discussed 
below. 


Path  (1)  satisfies  both  differential  expressions  (1)  and  (2)  at  the  expense  of 
(14).  Finally,  path  (2)  enforces  the  v-definition  (2)  and  the  difference 
system  (14)  at  the  expense  of  the  ODE  system  (1). 

Remark.  It  should  not  be  surmised  that  the  flow  chart  of  Figure  1 exhausts 
all  possible  computational  paths.  There  is  in  fact  another  path  family  in 
which  steps  (c,d)  are  carried  out  prior  to  (a,b).  The  practical  implemen- 
tation of  these  sequence  types  is  tied  up,  however,  to  very  specific  choices 
of  A and  B,  which  become  (linear)  functions  of  the  stepsize  h.  A study  of 
these  alternative  implementations,  which  offer  certain  organizational 
advantages,  will  appear  in  a sequel  article. 

Selection  of  Auxiliary  Vector 

The  behavior  of  the  propagated  computational  error  is  independent  of  the 
selection  of  auxiliary  vector  v through  Eq.  (2),  but  the  computational 
effort  per  step  is  not.  For  certain  choices  of  A and  B the  volume  of  com- 
putations--and,  to  a lesser  extent,  the  storage  requirements  - -can  be 
significantly  reduced  because  some  of  the  computational  steps  displayed 
in  Figure  1 can  be  either  simplified  or  bypassed  entirely.  The  actual 
operation  counts  turn  out  to  depend  on  sparseness  characteristics  of  the 
matrices  M and  D,  as  well  as  on  certain  properties  of  the  integration 
operators,  and  are  discussed  in  a later  section.  We  proceed  next  to 
examine  some  computationally  convenient  choices  of  A and  B,  which  will 
be  henceforth  identified  as  formulations  or  "forms"  for  short. 

Conventional  Formulation  (C-Form) 

The  choice  of  y described  in  standard  textbooks -- e.  g.  , [1,2,4,  5]--and 
in  most  papers  dealing  with  second-order  systems  is 

y = u (19) 

which  corresponds  to 
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(20) 


= M , B = C = 0 (20) 

The  temporal  derivative  v is  simply  the  acceleration  vector  u.  With  this 
choice,  the  components  of  the  algebraic  system  (15)  become 

E = M + h^  D + h.  h.  K 

<-w  o ^ 0 p 

(21) 

g_  = (M  + h.  D)  h'^  + h.  M h'^  + h.  h-  f 
^ n n 6 p-n 

The  sequence  of  calculations  corresponding  to  paths  (O'),  (1)  and  (2)  of 
Figure  1 are  listed  in  Table  1,  where  they  are  identified  as  (CO'),  (Cl) 
and  (C2),  respectively. 

The  fact  that  the  choice  (19)  effectively  compresses  the  four  state  vectors 

(u  . u , V , V ) into  three  (u  , u , u ) leads  to  some  simplifications  in 
'^-n 

the  general  computational  sequences  of  Figure  1.  Some  of  the  steps  become 
trivial,  but  have  been  left  in  Table  1 for  compatibility  with  the  arrangement 
of  Figure  1.  Some  practical  considerations  pertaining  to  the  computational 
sequences  of  Table  1 follow. 

1.  (C  O')  Form.  The  determination  of  accelerations  in  step  (a) 
requires  that  the  mass  matrix  M be  nonsingular.  This  is  an 
unfortunate  restriction,  as  many  commonly  used  mass -lumping 
procedures  in  finite  element  programs  assign  zero  mass  to 
non -translational  degrees  of  freedom  (the  assignation  of  small 
"fictitious"  masses  does  not  solve  the  problem  because  M remains 
ill-conditioned).  If  M is  nondiagonal,  step  (a)  is  computationally 
exp  ensive. 

2.  (Cl)  Form.  The  appearance  of  M‘^  is  avoided  by  working  with  the 
inertial  force  term  M vi  throughout;  consequently,  M can  be  singular. 

3.  (C2)  Form.  The  requirement  of  nonsingular  M applies  only  to  the 
starting  procedure,  during  which  the  initial  acceleration 


V = u = M'^  (f  - D il  - K u ) 
— o — O ~ o ~ — O — o 


Li 
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Table  1 


COMPUTATIONAL  SEQUENCES  ASSOCIATED  WITH  THE 
CONVENTIONAL  (C)  FORMULATION  v = u 


a 

a„.l  = 

b,b' 

c 

h = h d - c 

— n -^1  -^1 

d 

1.^  - u ^ 

h - h b - a 

-n  -n  — n 

(CO') 

e 

• 

e ^M  + h D)  h“  + h„  M h f 

6~-n  P p 6-n 

f 

E = M + h,  D + h K 

1-^  /-w  0 P 0 

g 

h 

a 

Mu  ,-Du  ,-Ku 

^ -n-l  -Ti-1  — n-1  — n-1 

(Cl) 

b 

c 

u 1 = u , (trivial) 

— n-1  — n-1 

. . , u , jM  u . . u 

Mh  =hd  -Me 

~ — n 

d-h 

Same  as  (CO*) 

a 

u - = li  - (trivial) 
— n-1  —n-1 

b 

Skipped 

(C2) 

C-g 

Same  as  (CO*) 

A 

h 

• 

iL  = (iin  -iii;  -h„h);)/(h,  h.) 

— n p —n  p 0 

h 

• 

• - uU  , i_  •• 

“ b_  + h Q a 
— n p ^1 

5 


must  be  calculated  from  the  initial  conditions  u > u • The 

— o — o 

calculation  sequence  shown  in  Table  1 differs  slightly  from  that 
of  Figure  1 in  that  the  accelerations _u^  are  determined  (through 
a double-differencing  operation)  prior  to  the  calculation  of  velo- 
cities. The  sequence  in  Table  1 is  used  more  extensively  in 
practice  (especially  in  undamped  systems)  than  that  of  Figure  1 
but  the  numerical  behavior  of  both  implementations  is  identical. 

Note  that  formulation  (CO)  is  not  shown.  The  implementation  of  path  (0) 

always  demands  that  v and  u be  kept  as  separate  arrays  and  is  consequently 

incompatible  with  the  compression  to  three  state  vectors  (u  , u , vi  ) 

^ —n  — n ""n 

induced  by  the  choice  of  (19). 

Jensen's  Formulation  (J -Form) 

Jensen  [3]  has  advocated  the  selection 

A=I,  B = C = D (23) 


which  gives 


V = M u + D u 


v = Mii  + Du  = f-  Ku 

The  auxiliary  vector  v can  be  interpreted  as  an  array  of  generalized 
momenta.  Its  temporal  derivative  v represents  the  unbalanced  dynamic 
force,  sometimes  called  the  "impulse"  in  classical  mechanics. 

With  the  choice  (23),  Eqs.  (16)  become 


E = M + h^  D + h K 

^ ^ P U p 


M + h,  h f 


*n  ^ — n p—n  6 P“n 

If  the  damping  matrix  D does  not  vanish,  (25b)  is  slightly  simpler  than 
(21b). 
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b 

^-1 

= 

V • 

h , + h . V 
— n-1  6— n-1 

c 



, ,V  V 

h d - c 

— n 

— n 

d 

_ 

a-  u 

h b - a 

—n 

— n — n 

e 

= 

M h'^  + h„  f 

^ ‘-n  P — n 6 P-n 

f 

E 

_ 

M + D + h K 

^ P ^ 6 P 

^-1 

g 

u 

— n 

h 

• 

u 

— n. 

= 

K - 

Mu  • + D u , 
^ — n-1  — n -1 


(v  , - h J/h 
n-1  -n-r  6 


^-1 


Same  as  (JO) 


(J2) 


Table  2 


COMPUTATIONAL  SEQUENCES  ASSOCIATED  WITH  THE  JENSEN  (J) 

FORMULATION  v = M u + D u 


The  sequence  of  computations  corresponding  to  paths  (0)  and  (2)  of 
Figure  1 are  displayed  in  Table  2,  where  they  are  identified  as  the  (JO)  and 
(J2)  formulations,  respectively.  Formulations  (JO')  and  (Jl)  are  not  shown 
as  they  lack  any  practical  value  (the  verification  of  this  assertion  is  left 
as  an  exercise  to  the  interested  reader). 

In  both  J-forms,  the  presence  of  a singular  M does  not  cause  any  particular 
difficulty.  There  are  also  no  special  problems  associated  with  the  starting 
procedure. 
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In  order  to  make  an  a priori  comparison  of  the  computational  effort 
involved  in  specific  implementations,  it  is  necessary  to  make  some 
simplifying  assumptions: 


1.  The  computational  effort  per  step  required  to  integrate  a linear 

dynamic  system  with  fixed  h is  considered  to  be  dominated  by 
matrix/ vector  operations.  Operations  involving  only  vectors 
are  neglected.  Diagonal  matrices  are  considered  equivalent  to 
vectors.  The  assembly  and  factorization  of  E--step  (f)  in 
Figure  l--is  not  included  in  the  count,  as  this  operation  is  per- 
formed only  once  for  fixed  h and  its  effect  on  all  formulations  is 
identical.  Similarly,  the  effort  spent  in  starting  and  stepsize- 
change  procedures  is  ignored. 

Z.  The  unit  of  work  (U)  is  a stiffness -matrix/ vector  multiply,  e.  g.  , 

K u . Matrices  M and  D are  assumed  to  be  either  diagonal 

— n 

(D  may  be  null)  or  have  the  sparseness  structure  of  K.  It  follows 
that  the  E matrix  (Zla)  or  (Z5a)  has  the  sparseness  structure  of  K. 
The  solution  of  the  linear  system  (15)  is  counted  as  ZU  (one  unit 
for  the  forward  pass  and  one  for  the  backsolve).  This  value  is 
accurate  only  for  fairly  dense  systems,  but  is  not  far  off  for 
matrices  stored  in  a ’’skyline*'  or  variable -bandwidth  format. 

3.  In  the  case  of  large  discrete  systems  that  require  out-of-core 

processing,  the  number  of  work  units  gives  an  idea  of  the  number 
of  passes  of  block-partitioned  matrices  through  core  buffers. 

On  many  computer  systems,  these  I/O  activities  may  dominate  the 
total  run  cost. 


r 


4.  The  effort  involved  in  the  evaluation  of  the  applied  force  vector 

_f(t)  is  not  considered,  as  this  is  not  implementation  dependent. 

Table  3 summarizes  the  work  per  step  associated  with  the  implementations 
described  in  the  preceding  section.  Work  units  are  tabulated  as  a function 
of  commonly  occurring  sparseness  attributes  of  M and  D as  well  as  of 
certain  properties  of  the  integrator  discussed  in  the  following  subsection. 
Optimal  formulations,  in  the  sense  noted  in  the  table,  are  given  for  each 
case.  It  is  clear  that  the  (JO)  formulation  is  preferable  in  the  majority 
of  situations. 

Backward  Difference  Operators 

A backward  difference  (BD)  operator  is  one  possessing  no  "historical" 
temporal  derivative  terms.  For  instance,  p.  = 0 (i  = 1,  . , , m)  in  Eq,  (7a). 
BD  operators  represent  a subset  of  the  general  class  of  LMS  operators 
that  is  widely  used  to  treat  stiff  differential  systems.  Important  members 
of  that  subset  include  Gear*s  stiffly -stable  family  [4]  (which  is  A-stable 
for  m = 1,2),  and  the  3-step,  A-stabie  Park's  method  [9]. 

As  indicated  in  Table  3,  considerable  reductions  in  the  computational 
effort  per  step  can  be  realized  in  most  formulations  if  a pair  of  BD 
operators  is  used  in  Eq.  (7).  For  example,  in  the  (Cl)  and  (C2)  forms, 
the  calculation  of  acceleration  terms  is  no  longer  required.  The  bene- 
ficial effect  on  vector  storage  requirements  is  discussed  in  Appendix  B, 

Error  Propagation  Behavior 

The  following  considerations  represent  a synopsis  of  part  of  the  material 
expounded  in  Part  II,  which  is  selected  insofar  as  it  impacts  the  rating  of 
the  formulations  listed  in  Table  3, 

It  is  well  known  that  the  response  u(t)  of  the  linear  differential  system 
(1)  can  be  expressed  as  the  superposition  of  normal  components  pertaining 
to  the  frequencies  tu.  and  mode  shapes  x.  of  the  associated  vibration 
eigenproblem 
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Table  3 

COMPARISON  OF  COMPUTATIONAL  EFFORT  REQUIRED  BY  VARIOUS  FORMULATIONS 
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0 


(26) 


L 


((ju?M+u).  D + K\  X.  - 
\ \ 1 ^ — I 

A similar  decomposition  applies  to  the  local  computational  error  sources, 
and  to  the  propagation  and  accumulation  of  such  sources. 

Discrete  element  dynamic  models  of  mechanical  systems  generally  fall 

into  the  category  of  stiff-oscillatory ' differential  systems.  This  term  is 

used  to  mean  that  the  eigenfrequencies  lie  close  to  (or  on)  the  imaginary 

axis,  and  that  the  ratio  | I spans  several  orders  of  magni- 

A Q max  mm 

tude;  typically  lO"*  to  lo“.  If  these  systems  are  treated  by  implicit  time 
integration  procedures,  an  A-stable  scheme  [4]  is  normally  required. 

For  most  problems  in  linear  structural  dynamics,  the  energy  content 
characteristics  of  the  excitation ^(t)  - -i.  e.  , its  energy  frequency  spectrum-- 
is  such  that  the  response  u(t)  can  be  adequately  represented  by  the  super- 
position of  a set  of  normal  responses  pertaining  to  low  frequency  modes. 
Accordingly,  the  stepsize  h is  chosen^  so  that  the  circular  sampling 
frequencies  Q.  = h vary  from  Q « 1 (corresponding  to  normal  compon- 
ents to  be  accurately  traced)  to  Q » 1 (corresponding  to  normal  components 
to  be  "filtered  out"  by  the  integrator).  (A  more  precise  stepsize  classifica- 
tion scheme  is  discussed  in  Part  II.  ) 

The  propagation  of  a modal  error  source  to  downstream  solutions  depends 
on  three  factors:  the  sampling  frequency  Q,  the  computational  path 
followed,  and  certain  properties  of  the  integration  operator.  It  turns  out 
that  the  amplification  of  a local  error  can  only  be  significant  in  the  low- 
frequency  (Q  « I)  and  high-frequency  (0  » 1)  regions  of  the  modal 
spectrum. 

In  the  low-frequency  end  (Cl  « 1),  the  error  propagation  is  primarily 
affected  by  the  computational  path.  Pertinent  results  are  summarized  in 
Table  3.  If  path  (O)--or  its  variant  (O') --is  followed,  the  downstream 

The  qualifier  "stiff"  is  used  here  to  denote  a mathematical  property  (widely 
different  time  constants)  and  is  unrelated  to  mechanical  stiffness. 

^In  practice,  the  stepsize  h is  chosen  adaptively  by  the  integration  program, 
which  "senses"  the  characteristics  of  the  excitation  through  accuracy- 
monitoring procedures. 
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amplification  of  a modal  error  source  is  asymptotically  independent  of 
the  stepsize.  On  the  other  hand,  if  paths  (1)  or  (2)  are  followed,  local 
modal  errors  are  amplified  by  a factor  asymptotically  proportional  to 
0 . For  example,  if  Q = ojh  = 0.  01  (a  common  sampling  rate  for  the 

fundamental  frequency),  the  error  amplification  difference  between 
paths  (0)  and  (1,  2)  is  two  orders  of  magnitude,  and  this  result  holds 
regardless  of  the  particular  integrators  selected  in  (7). 

The  error  amplification  behavior  in  the  high-frequency  region  O » 1 is 
controlled  by  the  computational  path  followed  and  by  the  coefficients 
and  6^  of  the  integration  operators  (7).  As  high-frequency  error  propaga- 
tion is  relatively  unimportant  in  linear  structural  dynamics,  it  has  not 
been  included  as  a rating  factor  in  Table  3. 
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Section  6 

STARTING  PROCEDURES 


The  actual  implementation  of  the  starting  procedure  for  LMS  operators 
with  m ^ 2 is  rarely  mentioned  in  the  literature.  The  most  common 
procedure  consists  of  using  a one-step  method,  such  as  the  trapezoidal 
rule,  to  generate  (m-1)  initial  solutions;  after  that  the  program  switches 
to  the  standard  integrator.  This  approach  usually  requires  one  extra 
factorization,  because  the  coefficient  matrix  (I6a)  of  the  starter  does 
not  necessarily  coincide  with  that  pertaining  to  the  integrator  (7).  We 
describe  here  a technique  that  avoids  the  additional  factorization. 


Introduce  the  ’’starting  family"  of  L.MS  operators: 


+ 

u 

a. 

= h_ 

+ h b^ 

— 1 

-1 

p 

-1 

—1 

u 

a^ 

= h 

. 

■t  h b^ 

—2 

-2 

P 

—2 

-2 

(1-step) 

(2-step) 


(27) 


u +a  =h^u  +hb 
— m — m P — m — m 


(target  operator) 


which  are  characterized  by  identical  values  of  the  generalized  timestep 

hp.  A similar  family  is  constructed  for  ^ if  h^  / h^  . As  an  example, 

if  the  target  integrator  for  u is  the  3-step  Park’ s method  [9]  , for  which 

h = 0.  6 h,  a convenient  starting  family  is 
P 


- 

u 

— o 

~ 

0.  6 h 

* 

0.  4 h u 
— o 

1.  2 

^1 

+ 

0. 

= 0.  6 h u^  t 

u - 
— n 

1.  5 

-1 

-H 

0,  6 u „ - 0.  1 u 

— n-2  — n 

-1 


(28) 


— n 


The  two  auxiliary  formulas  (28a,  b)  are  A-stable  [(28b)  is  obtained  by  com- 
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bining  40%  of  the  trapezoidal  rule  with  60%  of  Gear's  two-step  integrator]; 
hence*  no  spurious  high-frequency  oscillations  can  be  triggered  by  the 
starter. 

A change  in  the  stepsize  h requires  an  additional  assembly  and  factoriza- 
tion of  E in  most  existing  implementations*  A matrix -scaling  technique 
can  be  used,  however,  to  eliminate  such  refactorizations  in  linear  or 
nonlinear  problems  [10]  . 


f 

\ 


] 

I 

i 

i 

j 


1 


-1 


Section  7 
CONCLUSIONS 

A general  framework  for  categorizing  computer  implementations  of  direct 
time  integration  procedures  based  on  the  treatment  of  an  equivalent  first- 
order  system  by  implicit  LMS  operators  has  been  presented*  Five  parti- 
cular Implementations  associated  with  computationally  convenient  choices 
of  auxiliary  vector  and  with  certain  calculation  sequences  have  been 
studied  in  some  detail.  It  is  appropriate  to  summarize  in  this  final  section 
the  relative  merits  and  disadvantages  of  those  specific  formulations. 

The  (CO')  form  has  excellent  error  propagation  characteristics  for  low-  ; 

frequency  components,  but  is  substantially  more  expensive  than  the  other  j 

forms  if  either  M or  D is  nondiagonal.  The  explicit  determination  of 
accelerations  requires  the  presence  of  a well-conditioned  mass  matrix*  * 

As  this  form  is  in  no  case  superior  to  (JO),  it  is  not  recommended  for 
practical  use. 

Both  (Cl)  and  (C2)  display  poor  error  propagation  behavior  as  regards  low- 

frequency  response  components,  since  local  errors  are  magnified  by  an 

amount  inversely  proportional  to  the  stepsize  h*  (It  should  be  noted  that 

for  many  of  the  integration  operators  studied  in  Part  II,  this  amount  is 

larger  if  the  (C2)  form  is  used. ) The  (Cl)  form  has  the  advantage  of  not 

being  affected  by  a singular  (or  ill-conditioned)  mass  matrix,  because  \ 

the  product  M u can  be  kept  throughout  the  advancing  cycle.  Implementa-  i 

I tions  based  on  (C2)  may  run  into  starting  problems  if  M is  singular  or  ! 

1.  ^ ? 
[ ill-conditioned,  as  is  often  the  case  in  finite  element  models  that  include  < 

t j 

, rotational  degrees  of  freedom.  If  a BD  operator  pair  is  used  in  (7),  these  | 

two  forms  are  computationally  equivalent,  and  the  preceding  remarks  do 

! not  apply. 

1 

I 

The  (JO)  form  emerges  as  an  easy  winner  for  a general  purpose  implemen- 

I 

tation.  It  has  excellent  error  control,  is  not  affected  by  the  condition  of 


1-26 

L 


the  mass  matrix,  and  minimizes  the  computational  effort  per  step  in 
the  general  damped  case. 

The  (J2)  form  can  be  considered  a convenient  adjunct  of  (JO),  It  is  easy 
to  implement  both  (J)  forms  in  the  same  computer  program,  since  their 
storage  arrangements  (cf.  Appendix  B)  are  similar,  and  a simple  branch 
takes  care  of  the  difference  in  steps  (a,  b)  of  Table  2.  The  (J2)  form  is 
computationally  optimal  in  the  frequent  case  of  an  undamped,  lumped- 
mass  model,  and  could  be  used  in  such  instance  if  the  analyst  is  not  parti- 
cularly concerned  about  error  propagation  characteristics. 


] 


] 

i 

i 
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Appendix  A 

TWO -DERIVATIVE  METHODS 


In  the  discussion  of  numerical  integration  procedures  for  second-  and 
higher -order  ODE  systems,  it  has  been  customary --see,  e.  g, , [lf4]-- 
to  categorize  LMS  methods  as  follows: 

(I)  The  original  equations  are  recast  into  a first-order  system,  which 
is  treated  by  standard  (one-derivative)  DMS  operators.  This  is 
the  approach  followed  in  the  main  body  of  the  paper. 

(II)  DMS  formulas  containing  higher -order  derivatives  are  applied 
directly  to  the  higher-order  system. 

For  the  second-order  dynamic  equations  (1),  approach  (II)  requires  the 
introduction  of  an  L.MS  operator  containing  the  accelerations  u.  A 
general  form  of  such  schemes  has  been  studied  by  Geradin  [11]  . His 
expression,  rewritten  slightly  to  conform  to  the  notation  style  of  Eq.  (7), 
is 


where  7]  is  a scalar  used  to  effect  the  normalization  = I in  (Alb),  and 
^ 0 for  implicitness.  The  identification  of  the  3.  coefficients  in  both 
right-hand  sides  follows  from  elementary  consistency  considerations. 

Explicit  expressions  of  the  Houbolt,  Newmark,  and  Wilson  operators  in 
the  multistep  form  (Al)  may  be  found  in  [11]  ; some  specific  examples  are 
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given  in  Appendix  C (after  Part  II).  Eqs.  (Al)  may  be  rewritten  using  a 


compact  notation  similar  to  Eqs.  (9): 


u 

— n 

+ 

u 

a 

— n 

h 

—n 

u 

— n 

+ 

u 

c 

— n 

+ Tlh^ 

• • 

b'^ 
— n 

(A  2) 


Substracting  (A2a)  from  (A2b)  to  eliminate  and  collecting  historical 
terms,  the  LMS  pair  (A2)  may  be  presented  in  a form  analogous  to  (12): 


V 

— n 


h d \ 

p-r 


h-  V 
6 — r 


> + 


h'' 

— n 


(A3) 


where 


= u 
— n 


V 

— n 


u 

-n 


— n 6 ' 


(A4) 


Substitution  of  (A3)  into  (16)  gives 


E = 

/-w 

M 

+ h D 
0 ~ 

+ *'b'’6 

K 

M , 

(h.  h b'^ 
1 6 -n 

- + 

D 

+ h„  h f 
P 6-n 

(A5) 


If  D = 0,  the  algebraic  system  (15)  reduces  to 

IM  + “b ''e  5>  H..  “ ^ '‘e'’64 


(A6) 


which  expresses  the  well-known  fact  that  the  approximations  of  velocities 
through  (Ala)  is  not  required  for  the  direct  integration  of  an  undamped 
system. 

The  implementation  of  these  methods  is  seen  to  be  inextricably  linked  to 
the  conventional  choice  (19)  of  auxiliary  vector.  Accordingly,  the  three 
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computational  paths  (O',  1,  2)  of  Table  1 may  be  used  (with  appropriate 
modifications  to  account  for  the  expression  (A5b)  of  ^).  But  it  is  not 
difficult  to  show  that  (CO')  and  (Cl)  are  computationally  identical  in  the 
case  of  an  undamped  system,  and  very  similar  if  there  is  slight  damping. 
Hence,  only  two  implementation  variants  need  to  be  considered:  (CO' /Cl) 
and  (C2). 

Both  implementations  turn  out  to  display  local  error  amplification  of  order 
h ^ for  low-frequency  components,  i.  e.  , the  beneficial  effect  of  the  (CO') 
form  on  error  propagation  control  is  lost  if  a two -derivative  operator  (Alb) 
is  used.  As  (CO' /Cl)  consistently  involves  equal  or  larger  computational 
effort  per  step  than  (C2)--cf.  Table  3--the  latter  emerges  as  a clear  choice. 

Henrici  has  described  [1,  Ch.  6]a  radical  modification  of  (Alb),  called  the 
"summed  form,"  which  eliminates  the  0(h  ^)  dependence  of  the  propagated 
local  error.  For  an  undamped  (D  = 0)  system  (1),  this  device  is  equivalent 
to  using  the  reduced  first-order  form: 


(Mu') 

[ M V ) 

1 1 

= 

Eq.  (A7a)  is  treated  by  a standard  LMS  operator  derived  from  (Alb)  by 
formal  division  of  the  polynomial  Yq  + . . . + Y^  W (z-1)  while  (A7b) 
is  treated  by  a backward-Euler  operator.  This  technique  essentially 
amounts  to  a splitting  of  (Alb)  into  the  formal  product  of  two  one -derivative 
operators.  Therefore,  effective  control  of  computational  error  propagation 
demands  a reduction  to  a first-order  system,  whether  explicitly  done  as 
in  (3)  or  concealed  under  a special  name. 
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Appendix  B 

AERAY  STORAGE  REQUIREMENTS 


[ 


This  Appendix  should  be  of  interest  to  implementors  involved  in  the  pro- 
gramming of  one  of  the  formulations  defined  in  Tables  1 and  Z.  Table  4 
provides  a "snapshot"  view  of  array  structures  involved  in  forms  (Cl) 
and  (JO),  which  are  chosen  to  typify  the  (C)  and  (J)  formulations, 
respectively.  The  flow  of  computations  required  to  advance  the  solution 
over  the  current  (n-th)  step  is  indicated  with  arrowed  paths.  Matrix 
storage  is  not  explicitly  shown,  as  this  is  formulation -independent. 

If  the  number  of  degrees  of  freedom  (n£)  does  not  exceed  a few  thousands, 
all  vectors  displayed  in  Table  4 can  be  kept  in  high-speed  memory  (core) 
throughout.  Vector  storage  counts  in  such  case  are 


Form 

General  Operators 

BD  Operators 

(Cl) 

3 m + 3 

2 m + 1 

(JO) 

4 m + 4 

2 m + 2 

where  account  is  taken  of  the  fact  that  some  of  the  vectors  shown  in 
Table  4 can  be  overwritten  as  the  computations  proceed.  For  most  LMS 
operators  used  in  practice,  m = 1 to  3,  which  gives  vector  counts  ranging 
from  3 to  16.  The  beneficial  effect  of  using  BD  operators  is  clear. 

If  n£  exceeds  a certain  value  (typically  3000  to  10000),  it  is  necessary  to 
place  some  or  all  vectors  on  auxiliary  storage.  In  this  case  the  minimal 
number  of  core  vectors  becomes  2,  one  of  which  serves  as  a read  buffer. 

In  a dynamic  memory  management  environment,  the  transient  analyzer 
typically  begins  by  requesting  storage  for  matrix  block  buffers  to  the 
data  management  system.  Remaining  core  storage  is  then  parsed  in  n£- 
word  blocks  to  accommodate  vector  structures.  The  vector  storage  area 


i 


] 
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Table  4 

VECTOR  STORAGE  ARRANGEMENT  FOR  (Cl)  AND  (JO)  FORMULATIONS 


Multi  step  Historical 
Information  (missing  if  m=l) 


Previous  Step 
Information 


Current  Step 
Information 


u 

— n-m 


u 

—n-m 


[M  u 

—n-m 


u , 
— n-2 


Mu  J 
^ —n-2”* 


^ 1 
^ -n-1 


! ^-1 


^ — n-1 

f 


/f 


1.  If  a BD  operator  pair  (b  = d = 0)  is  used,  information  enclosed  in  brackets 

is  not  required  for  advancing  the  computations,  and  those  computational  steps 
identified  by  broken-line  paths  may  be  omitted. 

Z.  Storage  arrangement  for  M,  D,  K,  E is  not  shown,  as  matrix  storage 
requirements  are  formulation -independent. 
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] is  commonly  administered  as  a "revolving  pool,"  which  can  be  viewed  as 


a rectangular  array  of  row  dimension  n^  whose  columns  are  assigned  to 
specific  vectors  through  a pointer  string  linked  to  a cyclic  reservation 
scheme* 
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ABSTRACT 


The  propagation  of  computational  error  in  the  direct  time  integration 
of  the  equations  of  structural  dynamics  is  investigated.  Asymptotic 
error  propagation  equations  corresponding  to  the  computational  paths 
presented  in  Part  I are  derived  and  verified  by  means  of  numerical 
experiments.  It  is  shown  that  there  exists  an  implementation  form  that 
achieves  optimum  error  control  when  used  in  conjunction  with  one- 
derivative  methods.  No  such  form  is  found  for  two-derivative  methods. 
A numerical  beating  phenomenon  is  observed  for  certain  implementa- 
tions of  the  average  acceleration  method  and  the  trapezoidal  rule,  which, 
from  an  error  propagation  standpoint,  is  highly  undesirable. 
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Section  1 
INTRODUCTION 


The  material  presented  in  Part  I may  be  sufficient  for  those  readers  who 
are  primarily  interested  in  organizational  aspects  of  the  computer  imple- 
mentation of  time  integration  procedures.  Part  II,  on  the  other  hand,  is 
devoted  to  a detailed  analysis  of  the  computational  error  associated  with 
the  calculation  sequences  described  in  Part  I,  and  provides  substantiation 
for  claims  made  as  regards  the  error  propagation  behavior  of  those 
implementation  forms.  This  Part  endeavors  to  answer  questions  such  as: 
which  implementation  form  minimizes  error  propagation  effects?,  what 
are  the  advantages  of  one-derivative  LMS  operators  over  two-derivative 
methods?,  how  do  implementation  forms  affect  the  algorithmic  properties 
of  the  integrator?,  and,  finally,  how  can  the  importance  of  computational 
errors  be  assessed  in  advance  of  the  actual  computations?  The  answers 
to  these  questions  have  obvious  bearing  on  the  implementation  of  advanced 
error  control  strategies  in  time  integration  packages  developed  to  support 
general-purpose  structural  analyzers. 

The  potential  importance  of  computational  error  propagation  can  be 

appreciated  from  the  following  considerations.  A transient  response 

Z 4 

analysis  by  implicit  integration  procedures  typically  involves  10  to  10 
time  steps.  Each  step  in  turn  requires  a substantial  number  of  algebraic 
manipulations.  Computational  errors  introduced  at  each  time  step  are 
propagated  to  subsequent  time  stations  by  the  integration  operator.  The 
compound  effect  of  propagation  and  accumulation  of  local  errors  can 
seriously  jeopardize  the  accuracy  of  the  solution  in  many  instances, 
especially  in  large-scale  problems. 

Organization  of  Part  II 

A summary  of  new  results  obtained  in  this  Part  is  provided  to  aid  the 
reader  in  progressing  through  subsequent  sections.  An  introductory  sec- 
tion presents  basic  concepts  required  in  the  error  analysis.  Sources  of 


i 

j 
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local  computational  error  are  exhibited  in  some  detail,  and  mechanisms 
of  propagation  and  accumulation  described.  Governing  equations  for  the 
propagation  of  a single  local  error  source  for  a linear  system  are  then 
introduced;  it  is  shown  that  the  study  of  the  propagated  error  suffices  to 
characterize  the  behavior  of  the  overall  computational  error. 

Computational  error  sources  are  classified  into  operational  and  inherent 
according  to  the  functional  dependence  of  the  associated  propaged  error 
on  the  stepsize  h.  The  operational  error  is  introduced  by  inexactness  in 
the  solution  of  the  algebraic  equations  that  must  be  solved  at  each  step, 
whereas  the  inherent  error  is  due  to  uncertainties  in  the  specification  of 
applied  forces. 

A linear,  undamped  oscillator  model  is  used  to  examine  the  behavior  of 
the  propagated  operational  error  as  a function  of  the  sampling  rate.  An 
asymptotic  analysis  technique  is  used  to  show  that  the  long  term  behavior 
of  the  propagated  operational  error  can  be  characte*ri'  ^d,  for  small  step- 
sizes,  by  a limit  differential  or  integro-differential  propagation  equation. 

Examination  of  the  amplification  factor  associated  with  these  equations 
indicates  that  the  computational  sequence  labeled  as  path  (0)  in  Part  I 
is  immune  to  local  error  propagation,  hence  achieving  optimal  accumulated 
error  behavior. 

An  important  source  of  local  inherent  error  in  nonlinear  structural 

dynamics  is  exhibited  by  consideration  of  a sample  nonlinear  system  treated  , 

1 

by  the  pseudo-force  linearization  technique.  While  a complete  analysis  of  ^ 

the  inherent  error  propagation  in  nonlinear  problems  is  not  undertaken,  its  -j 

qualitative  "semilocal"  behavior  is  illustrated  through  a linearized  treat-  i 

ment.  A more  thorough  treatment  of  the  coupling  of  the  propagated  inherent 
error  with  computed  solutions  of  highly  nonlinear  problems  will  be  the  sub-  | 

ject  of  future  investigations. 

I 

I 

Numerical  experiments  to  characterize  theerror  propagation  behavior  of 
several  commonly  used  integrator  operators  have  been  performed  over  a j 

five -magnitude  range  of  sampling  frequencies  cuh,  and  the  results  are  ^ 

displayed  in  the  form  of  amplification  factor  spectra.  In  the  low-frequency  •, 

end,  these  results  agree  very  closely  with  the  predictions  of  the  asymptotic  i 
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analysis.  In  the  high-frequency  end,  an  interesting  error  growth  is  noted 
for  weakly  A-stable  operators  such  as  the  trapezoidal  rule  (and  its  two- 
derivative  counterpart,  the  average  acceleration  method)  under  certain 
computational  paths.  This  phenomenon  is  shown  analytically  to  be 
caused  by  "numerical  beating"  associated  with  the  asymptotic  appearance, 
as  u)h  -*  00 , of  secular  terms  related  to  closely  allied  characteristic 
roots  of  the  propagation  equation.  Finally,  the  implications  of  the  pre- 
ceding results  on  the  organization  of  problem -adaptive  time -integration 
software  are  discussed. 


Section  2 

SUMMARY  OF  NEW  RESULTS  IN  PART  II 


The  analysis  of  computational  error  propagation  in  linear  dynamic 
problems  relies  on  the  decomposition  of  error  sources  over  a frequency 
spectrum  encompassing  sampling  rates  normally  encountered  in  practice. 
To  our  knowledge,  this  approach  has  not  been  previously  used  for  that 
purpose.  More  specific  results  include; 


1.  Classification  of  computational  error  sources  according  to 
their  importance  with  respect  to  the  sampling  rate. 

2.  Derivation  of  deterministic  local  error  propagation  equations 
for  the  implementation  forms  discussed  in  Part  I. 

3.  Tabulation  of  accumulated  error  estimates  as  function  of 
implementation  form  and  local  error  correlation  properties. 

4.  Exhibition  of  a numerical  beating  phenomenon  associated  with 
high-frequency  error  propagation  of  certain  implementations  of 
the  trapezoidal  rule  and  the  average  acceleration  method. 

5.  Display  of  experimentally  generated  spectral  error  propagation 
characteristics  of  specific  integrators,  and  comparison  with 
analytical  predictions. 


} 
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Section  3 

COMPUTATIONAL  ERRORS 


Definition 


Let  u(t.)  denote  the  solution  of  the  difference  equations  (15)  obtainable 
with  exact  computations.  The  computational  error  _e(t.)  is  defined  as  the 
residual 

^ = £{L)  = 'u(t.)  - u(L)  (29) 

where  ^(t.)  is  the  computed  solution.  The  computational  error  results 
from  the  propagation  and  accumulation  of  local  errors  committed  at  each 
integration  step.  Consequently,  an  analysis  of  (Z9)  must  begin  with  a 
study  of  the  local  errors. 

Local  Errors 


Local  errors  may  be  studied  in  vacuo  by  ignoring,  for  the  moment,  the 

processes  of  propagation  and  accumulation.  At  each  time  step  t , we 

n 

effectively  solve  the  "perturbed"  equation 

(E  + aE)  + AUj^)  - ^ (30) 

instead  of 

= in  (‘5) 


In  Eq.  (30),  and  represent  the  errors  made  in  the  calculation  of 
the  right-hand  side  e and  solution  vector  u , respectively,  and  AE  is  an 
equivalent  perturbation  of  the  coefficient  matrix  E in  the  spirit  of  the 
backward  error  analysis  of  Wilkinson  [12]  . We  use  the  notation 


€ 


u 


^u 

-TX 


e 

g 


II  II 
II  ^11 


(31) 
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to  identify  relative  local  errors,  where 
vector  norm. 


denotes  any  appropriate 


Sources  of  Local  Error 


The  three  main  sources  of  local  error  mentioned  in  the  Introduction  to 
Part  I will  be  examined  now  in  more  detail.  For  the  ensuing  discussion 
we  refer  to  an  expanded  form  of  (30),  which  display  error  terms  pertain- 
ing to  difference  sources: 


E (On  + Au;  + ^ 


(32) 


in  which 


= -AE  + AU^)  -AE 


(33) 


and  where  other  undefined  terms  are  defined  in  the  text  below. 


Rounding  errors  caused  by  the  use  of  finite  precision  arithmetic  contami- 
nate all  of  the  computational  steps  exhibited  in  Tables  1 and  2.  The 
aggregate  effect  of  those  errors  is  collected  in  the  right-hand  side  per- 
turbation  , which  induces  a solution  error  Au^  . The  term  is 
relatively  unimportant  in  most  circumstances,  with  a possible  exception: 
errors  introduced  in  the  calculation  of  accelerations  in  step  (a)  of  the  (CO') 
form  (or  of  the  starter  of  the  (C2)  form)  if  the  mass  matrix  is  ill-conditioned. 

Solution  errors  made  in  solving  (15)  by  direct  or  iterative  techniques  are 

g 

represented  by  the  term  Au^  and  the  corresponding  right-hand  side  residual 
A^ . If  a direct  solution  technique  is  used,  the  matrix  perturbation  aE 
results  from  factorization  errors.  If  an  iterative  solution  procedure  is 

g 

used--as  required  in  nonlinear  problems--the  term  A^  is  augmented  by 
uncertainties  arising  from  the  acceptance  of  solution  iterates  that  have 
not  attained  the  full  digital  significance  allowed  by  the  condition  of  the 
iteration  operator. 
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r— 

Finally,  force  evaluation  errors  affect  the  accuracy  of  the  right-hand  side 
through  the  term  given  by  (33b),  in  which  denotes  the  uncertainty 

in  f . The  corresponding  solution  error  term  is  Au^  , 

—n  — n 

The  relative  magnitude  of  right-hand  side  and  solution  errors  is  related 
by  the  condition  number  of  matrix  E; 

e ^ C(E)  e 

u g 

(34) 

C{E)  = II  E|I  II  E"^I| 

From  the  error  analysis  viewpoint,  rounding  and  solution  error  terms  can 

be  lumped  together  into  a component  called  the  "operational  error"  source. 

2 

Because  of  the  presence  of  the  h factor  in  (3Zb),  force  errors  are  treated 
as  a separate  component  called  the  "inherent  error"  source. 

Remark.  The  term  "inherent  error"  was  introduced  by  Henrici  [1]  to 
identify  errors  in  the  evaluation  of  the  forcing  term  f(y,t)  while  integrating 
i the  single  ODE  y'  = f(y,t).  He  then  used  the  term  "induced  error"  to  de- 

^ note  the  equivalent  of  rounding  and  solution  errors;  we  prefer  the  more 

I suggestive  term  "operational  error." 

Error  Propagation 

The  foregoing  analysis  of  local  errors  would  be  sufficient  for  time- 
independent  (static)  problems.  In  a dynamic  system,  however,  local 
errors  committed  at  t.  are  propagated  to  the  stream  of  solutions  u(t,  ), 

j k 

i s j , by  the  feedback  effect  of  historic  terms.  The  overall  computational 
error  (28)  at  t^  may  be  expressed  as  the  accumulation  of  propagated 
\ sources: 


where  X.,  are  influence  coefficients  expressing  the  propagated  error  at 

J i 

t.  due  to  the  k-th  component  Au.,  of  the  local  error  Au.  at  t.  , and  G is 
I Jk  -j  j’  -n 

II-7 


a vector  operator  mapping  the  contribution  of  the  i-th  station  onto  e . 

i 

For  general  ODE  systems,  are  nonlinear  functions  of  the  solution 
vector  u.  In  the  case  of  a linear  system,  however,  the  principle  of 
superposition  can  be  invoked  to  eliminate  the  nonlinear  operator  , 
whence  (35)  simplifies  to 


X Au., 

Jk  Jk 


i=0  j=0 


where  the  are  independent  of  u,  and  can  be  collected,  after  summing 
on  i,  into  "propagation  matrices"  L.  For  constant  stepsize  h,  these 
matrices  assume  a regular  pattern  illustrated  below  for  the  case  n = 0,  ...  3: 


L,  1 

hz  h ^ 

^3  ^2 


where  L depends  only  on  the  station  index  difference  s = n - i.  It  is 
s 

therefore  sufficient  to  investigate  the  propagation  of  a typical  local  error 
source,  which  can  be  taken  to  be  for  conveniency.  The  resulting 
propagation  function  r(t  ) = L Au.  can  be  used  to  fill  up  the  first  column 
of  the  supermatrix  in  (37),  and  other  columns  follow  by  simple  shifts. 

The  accumulated  error  ^(t^)  can  then  be  calculated  by  applying  statistical 
assumptions  on  the  occurrence  of  local  errors  Au.  on  the  right-hand  side 
of  (37). 


A detailed,  component-by-component  calculation  of  e(t  ) is  seldom  useful 

— n 

(or  even  possible).  All  that  is  needed  for  practical  applications  is  a rough 
estimate  of  the  relative  error  (|  £l|/  ||  ^|[  given  the  local  level  of  accuracy 
(31)  and  some  gross  problem  identification  parameters.  An  important 
step  in  this  regard  is  to  effect  a spectral  decomposition  of  the  propagated 
error  by  passing  to  normal  coordinates,  and  to  identify  error  sources 
whose  frequency  spectrum  is  similar. 


k 


F 
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Operational  and  Inherent  Error 

For  reasons  outlined  above,  it  is  convenient  to  separate  the  propagated 
typical  local  source  into  two  components: 


L Au 
^n  — c 


= E(t„)  + 


(38) 


where  p(t)  denotes  the  operational  error  resulting  from  the  propagation  or 
rounding  and  solution  errors,  and  £^(t)  is  the  inherent  error,  which  results 
from  the  propagation  of  force  errors.  These  error  functions  satisfy  the 
linear  propagation  difference  equations 


n 


= 0,  1, 


(39) 


under  initial  conditions  discussed  in  following  sections.  The  right-hand 
side  terms  and  ^ are  to  be  evaluated  following  the  computational 
sequences  described  in  Part  I,  in  which  the  external  force  term  is 
set  to  zero. 


The  distinction  between  operational  and  inherent  error  arises  naturally  i 

j 

from  the  study  of  their  dependence  on  the  stepsize  h,  or,  more  precisely,  j 

the  circular  sampling  frequency  0.  = (jjh.  Broadly  speaking,  the  operational  | 

error  is  dominant  for  ’'small"  stepsizes  whereas  the  inherent  error  becomes  j 

important  for  "intermediate"  stepsizes.  The  concept  of  stepsize  magnitude 
is  defined  more  precisely  in  the  following  subsection. 

Modal  Decomposition 

As  noted  in  Part  I,  the  transient  response  of  a linear  dynamic  system  can 
be  viewed  as  the  superposition  of  normal  response  components  of  (generally 
complex)  frequencies  which  are  the  solutions  of  the  eigensystem  (26), 

Given  a specific  frequency  value  (U,  a stepsize  h will  be  called 

Small,  if  |u)h|  <0,1 
Intermediate,  if  0.1  <‘  | ojh  | < 2 
Large,  if  |cuh|  > 2 (Nyquist  frequency) 
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In  the  analysis  of  stiff-oscillatory  differential  systems  by  implicit  inte- 
gration procedures,  the  stepsize  h is  generally  such  that  there  are  normal 
components  in  the  three  foregoing  categories.  Generally  speaking,  small 
stepsizes  are  used  for  components  to  be  accurately  traced,  intermediate 
stepsizes  for  components  to  be  roughly  traced,  and  large  stepsizes  for 
components  to  be  "filtered  out"  by  the  integrator. 

All  error  vectors--local,  propagated  and  accumulated- -can  be  decomposed 
into  normal  components.  This  spectral  decomposition  technique  will  be 
used  in  following  sections  to  establish  practical  estimates  on  the  peak 
amplifications  of  the  propagated  operational  and  inherent  errors. 


n-10 


Section  4 

THE  OPERATIONAL  ERROR  PROPAGATION  PROCESS 


The  operational  error  is  dominant  for  those  modal  response  components 
associated  with  small  stepsizes.  Consequently,  this  is  the  only  error 
source  that  the  analyst  should  be  concerned  about  in  problems  where 
the  bulk  of  the  excitation  energy  is  contained  in  low  frequency  modes. 


As  noted  in  the  preceding  section,  it  is  sufficient  to  study  the  propagation 
of  a single  operational  error  source  at  t = 0,  and  the  resulting  propagation 
function  p(t)  can  be  projected  on  the  normal  coordinates; 


2.(t) 


ik 


(40) 


The  principal  features  of  a typical  modal  propagation  function  p(t)  = 
associated  with  a small  stepsize,  i.  e.  , \ u)^h|  < O.l,  are  illustrated  in 

Figure  2 for  the  undamped  case.  Three  stages  may  be  noted:  an  initial 
jump  due  to  local  error  initial  conditions,  a transient  (boundary  layer) 
period,  and  a free  oscillation  regime. 


Initial  Conditions 


The  appropriate  initial  conditions  for  the  modal  analog  of  (39a)  are 

P. 


* r . s 
= Au  + All 
o o o 


P “ P 

*^o  ^O  P 


P-l  = P-2 


(41) 


= 0 


where  the  value  of  p results  from  (12a),  in  which  h^  = 0 on  account  of 

o o 

(41c).  The  conditions  (41)  may  be  interpreted  as  the  application  of  a step 

'’error  velocity"  p at  t = -h  , as  depicted  in  Figure  2, 
o p 
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Transitional  Propagation  (TP)  Period 


The  TP  period  extends  from  t = 0 to  t = t„  , and  "bridges"  the  initial 

discontinuities  (40)  with  the  free  propagation  regime.  The  end  of  this 

period  may  be  conventionally  characterized  by  an  almost  uniform  velocity 

p = p(t  ) not  varying  by  more  than  1%  of  p , if  ujh  is  set  to  zero.  If  a 

one-step  method  is  used,  or  if  path  (0)  is  followed,  this  period  lasts  one 

step,  i.  e.  , t„  = h.  On  the  other  hand,  for  multiple  step  methods  (m  ^ 2) 
r 

implemented  using  paths  (1)  or  (2),  t„  depends  on  the  effect  of  the  initial 

r 

conditions  (41)  on  the  spurious  roots  of  the  integrator,  and  on  the  rapidity 

of  decay  of  such  spurious  oscillations;  nonetheless,  t_  2 m h represents 

r 

a good  approximation  for  most  commonly  used  A -stable  operators. 

The  exit  conditions  from  this  period  have  the  small-h  expansions 

= PPo^'^'^  + 0(1)  \ 

Pp  = ap^  + 0(Q)  / Q = o;h  (42) 

P„  = + o(n^)  / 

X o 

J.n 

p(t)  dt,  and  where  p , a and  Q are 
o 

scalars  of  order  unity.  The  values  of  p and  a , which  recur  in  the 
asymptotic  analysis  of  next  section,  are  tabulated  in  Appendix  C for 
several  commonly  used  integrators. 

Free  Propagation  (FP)  Regime 

For  t ^ tp  , the  propagated  modal  error  oscillates  with  a period  associated 
with  the  (generally  integrator -distorted)  frequency  pertaining  to  that  mode. 
The  amplitude  of  the  oscillation  is  primarily  determined  by  the  exit  con- 
ditions (42)  and  by  the  sampling  frequency  Q = ooh.  The  algorithmic 
properties  of  the  integrator  intervene  only  through  second-order  effects 
such  as  frequency  distortion  and  numerical  damping,  which  vanish  rapidly 
as  n approaches  zero. 

The  error  amplification  factor  is  defined  as  the  ratio 
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where  p is  the  largest  absolute  value  attained  by  p(t).  It  is  intuitively 

^max  V.  j 

clear  from  Figure  2 that  for  small  ujh  the  maximum  of  p(t)  is  reached  at 
the  first  amplitude  peak  in  the  FP  regime,  as  succeeding  maxima  are 
abated  through  numerical  damping  effects. 


Amplification  Analysis  Procedure 


The  central  problem  in  the  study  of  operational  error  propagation  consists 
of  establishing  the  dependence  of  the  error  amplification  factor  (43)  on  the 
sampling  frequency  for  specific  operators  and  implementation  procedures. 
A comprehensive  study  spanning  a wide  wh  range  requires  the  numerical 
integration  of  the  governing  difference  equation  (39a)  on  the  computer. 

This  is  in  fact  the  approach  followed  in  the  numerical  experiments  reported 
later.  Valuable  insight  on  the  effect  of  the  computational  path  on  (43)  for 
tuh  « 1 can  be  obtained,  however,  by  an  asymptotic  technique. 


For  a modal  component  of  frequency  to,  the  dependence  of  on  the 
stepsize  can  be  expressed  as  the  Laurent  series  in  0 = u)h: 


H 

P 


C 

o 


Q + . . . 


(44) 


where  c c ....  are  weighting  coefficients  depending  on  the  integrator 
and  computational  path  used.  The  purpose  of  the  asymptotic  analysis 
carried  out  in  the  following  section  is  to  determine  the  leading  (first  non- 
vanishing) term  in  (44).  The  analysis  involves  two  stag^^s. 

1.  The  calculation  of  the  TP  exit  conditions  (42)  by  solving  the 

propagation  difference  equations  (39)  up  to  tp  . The  leading 
expansion  terms  in  (42)  can  be  obtained  by  setting  h - 0 
ab  initio,  which  simplifies  the  computations  considerably. 

Results  for  specific  operators  ana  computational  paths  are 
collected  in  Appendix  C • 


The  asymptotic  behavior  of  p(t)  in  the  FP  region  (t  ^ tp)  is 
established  by  replacing  the  propagation  difference  equation  (39a) 


k 


by  the  differential  (or  integro-differential)  form  approached  as 
u)h  tends  to  zero.  A closed  form  solution  matching  the  TP  exit 
values  (4Z)  as  initial  conditions  at  tp  then  yields  the  desired 
information  on  peak  values. 

The  algebraic  manipulations  are  simplified  by  ignoring  the  damping  term 
D ^consistently  throughout  the  error  analysis,  and  working  with  the 
homogeneous  system 

Mp  + Kp  = 0 (45) 

or  its  modal  counterpart 

p + p = 0 (46) 

as  source  differential  systems  for  the  propagation  difference  equation 
(39a).  The  introduction  of  light  modal  damping  in  (46)  in  fact  trims  the 
amplification  factor  by  an  amount  of  the  order  of  the  first  neglected  term 
in  (44). 
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Section  5 

ASYMPTOTIC  ANALYSIS  OF  OPERATIONAL  ERROR  PROPAGATION 


The  variation  of  the  propagated  operational  error  in  the  FP  regime  can  be 
studied  as  an  undamped,  free-vibration  initial-value  problem.  In  this 
section  we  examine  the  limit  differential  forms  associated  with  computa- 
tional paths  (Z)  and  (0),  as  typified  by  the  (CZ)  and(J0)  formulations, 
respectively. 


(CZ)  Form 


The  error  propagation  equations  are  obtained  by  specializing  the  calcula 
tions  of  Table  1 to  Eq.  (45): 


(M+h.h  K)p  = M (h^  + h h^) 


or 


M 


p - h^  - hp  h^ 

^ --n- — + kp  = m 


^ + KPn  " ° 


(47) 


(48) 


where  the  vector  expression  postmultiplying  M has  been  recognized  to  be 
the  difference  approximation  used  for  the  accelerations  in  the  (CZ)  form. 
Consequently,  as  u;h  - 0,  (48)  approaches  the  differential  equation  (45), 

Its  general  solution  is  easily  expressed  by  passing  to  normal  coordinates. 
The  solution  of  the  associated  modal  equation  (46)  with  the  initial  conditions 
p^,  Pp  at  t = t^.  is 

p(t)  = Pp  cos  uu  (t  - tp)  + a;  ^ Pp  sin  a;  (t  - tp)  (49) 

Substituting  the  expansions  (4Za,  b)  into  (49)  and  estimating  the  peak  value 
gives 

H = p (u;h)'^  + 0(1)  (50) 
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(JO)  Form 


The  analysis  of  this  form  is  more  involved.  The  error  propagation  equa- 
tions follow  by  specializing  (25)  to  the  undamped  case: 


M hP  + h 
-n  6 — n 


(51) 


in  which  v = M p denotes  here  the  auxiliary  vector  (24a)  associated  with 
£(t).  This  vector  is  determined  recursively  from  the  differential  equation, 

i,  e.  , 


v 

— n 


- K 


En 


= h'' 
— n 


h^  v 
6 “n 


(52) 


To  eliminate  v from  (51),  the  historical  vector  ^ has  to  be  expressed 
in  terms  of  the  value  of  at  past  stations.  It  can  be  shown  that 


(53) 


in  which  pi.  are  coefficients  resulting  from  the  Taylor  expansion  of  the 
generating  rational  function 


6.  + 6 z + ...  6^  z 

o V m 


m 


V + V,  z + . . . Y z 
^o  U ^m 


m 


- ^Ao  + M'l  + • • • + 2 


n-1 


(54) 


where  ^.re  coefficients  of  the  ^-integrator  (7b).  The  -oo  limit  in 

(53)  can  be  replaced  by  zero  on  account  of  (41c),  Define  now 


En 


K‘^  h'" 


(55) 


It  turns  out  that  P is  a discrete  integral  of  £_(t).  This  can  be  verified  by 
expressing  (9)  for  and 


= h 


6En 


+ h 

— n 


cP 

— n 


(56) 


11-17 


in  which  and  ^ can  be  recursively  calculated  in  terms  of  using 
(52).  As  an  example,  for  the  trapezoidal  rule  (y^  = -1,  ~ ^1  " 


obtain 


= 1/2  h(£^  +2£^  + ....  +2£^_^  + 


which  is  the  well-known  trapezoidal  quadrature  formula.  Inserting  (53) 
into  (51)  produces  the  difference  expression 

M (p  - h^)/h  + K i = 0 

( 

or 

M p ^ ^ 

^ i-n  — 

where  ^ is  the  difference  expression  used  in  step  (b)  of  Table  2.  In  the 
limit  u)h  -♦  0 , (58)  becomes  the  integro-differential  equation 

t 

n 

M^+K^^dt^O^  ( 


The  solution  of  the  associated  normal  equation 

t 


p + uj  / pdt  = 0 


satisfying  the  initial  conditions  p^^  , Pp  at  t^  , is 

p(t)  = cos  uo  (t  - t^)  - oj  p^  sin  (w  (t  - t^)  (61) 

Introducing  (42b,  c)  into  (61)  and  estimating  the  peak  value  produces 

K = a + O(ujh)  (62) 

P 

The  error  amplification  factor  of  the  (JO)  form  is  asymptotically  independ- 
ent of  the  stepsize  h.  Consequently,  this  implementation  is  immune  to 
local  propagation  error  effects. 


Remark  _1.  The  choice  of  the  (J)  and  (C)  formulations  to  exhibit  the  error 
propagation  behavior  of  the  (0)  and  (2)  paths,  respectively,  has  no  special 
significance.  In  fact,  the  same  analysis  (leading  to  identical  conclusions), 
could  be  carried  out  for  an  arbitrary  auxiliary  vector. 
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Remark  2,  An  analysis  of  computational  paths  (O')  and  (1)  leads  to  the 
limit  forms  (59)  and  (45),  respectively.  The  corresponding  amplification 
factors  have  expansions  similar  to  (62)  and  (50),  respectively,  but  the 
coefficients  p and  a are  generally  different  (cf.  Appendix  C).  Conse- 
quently, a detailed  analysis  of  these  paths  will  not  be  given. 

Remark  3.  The  integro -differential  equation  (59)  that  governs  the  error 
propagation  behavior  of  the  (0)  path  has  the  important  property  of  being 
insensitive  to  the  initial  step  velocity  (41b).  This  smoothing  property 
has  been  actively  exploited  in  the  field  of  discrete  sequential  estimation 
for  the  design  of  stochastic  filters  [13].  The  elimination  of  error  ampli- 
fication effects  is  in  fact  analogous  to  the  filtering  of  process  noise 
through  a Kalman  filter, 

Two-Derivative  Methods 


An  analysis  of  the  operational  error  propagation  for  the  two -derivative 
operators  mentioned  in  Appendix  A leads  to  modal  amplification  factors  of 
the  form  (50),  regardless  of  the  computational  path  used.  The  smoothing 
effects  of  the  auxiliary  vector  sequence  (52)  are  lost  in  the  (0)  path 
because  computed  solutions  are  independent  of  the  velocities  in  the 
undamped  case.  Consequently,  the  error  filtering  effect  of  the  (CO)  and 
(JO)  implementations  is  lost  if  a two-derivative  operator  is  used.  This 
can  be  a significant  drawback  of  two -derivative  methods  when  the  opera- 
tional error  can  severely  affect  the  quality  of  the  answers. 
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Section  6 

ACCUMULATED  OPERATIONAL  ERROR 


The  accumulated  operational  error,  denoted  results  from  the 

^ f P 

composition  of  propagated  operational  sources.  An  estimation  of 
necessarily  requires  some  assumptions  on  the  statistical  distribution  of 
sequential  local  error  occurrences.  The  reader  is  referred  to  Henrici 
[1,  2]  for  a detailed  exposition  of  statistical  analysis  techniques.  We  are 
here  content  to  illustrate  the  accumulation  process  pertaining  to  a small- 
stepsize  modal  component  (uoh  < 0.  1)  in  the  simplest  instance  where  local 
modal  errors  are  perfectly  correlated,  i.  e.  , 


Au 


All. 

i 


for  each  i 


(63) 


For  computational  paths  (1)  and  (2),  the  dominant  term  in  the  propagation 
function  (49)  is,  for  small  Q = uoh  : 

p(t)  fts  u)"*  p^  sin  ua(t  - L) 


Consequ  ently, 


e 


P 


n 


-L  . 

Pp 


n 


sin  uo  (t 

n 


n 


i: 

;-n 


sin  jC 


= uj'  Pp  sin  l/2(n-l)f^  sin  1/2  0 cosec  1/2  Q (64) 

Expressing  p^  in  terms  of  (63)  through  (42a),  the  following  expansion  for 
the  peak  value  of  is  obtained: 

e^  = max  = p fl  ^ Au  + 0(Q  ) (65) 

max  n ‘ n‘ 


A similar  argument  for  path  (0),  using  the  dominant  term  in  (61),  yields 


e^  a Q ^ An  + 0(1) 

max 


(66) 


II -20 


Leading  terms  of  the  asymptotic  expansions  of  the  propagated  and  accumu- 
lated errors  in  the  low  and  high  frequency  regions  of  the  modal  spectrum 
are  collected  in  Table  5.  The  error  exponent  u defined  there  is  a measure 
of  the  correlation  of  sequential  local  error  occurrences. 


■Remark.  In  single  ODE's,  where  rounding  is  the  only  significant  local 

error  source,  it  is  generally  assumed  that  sequential  local  errors  are 

normally  distributed  about  a zero  mean,  which  leads  to  a favorable 

exponent  1/2.  For  large-scale  stiff  systems,  however,  that  assumption 

may  be  unrealistic.  If  the  system  is  at  least  moderately  ill-conditioned, 

the  solution  error  term  u tends  to  dominate  over  the  rounding  error. 

— n 

As  the  same  E is  normally  used  for  many  steps,  the  factorization  error 
perturbation  AL  remains  unchanged;  furthermore,  the  low-frequency  modal 
components  in  u^  vary  slowly  at  the  usual  sampling  rates.  The  net  effect 
is  that  successive  local  errors  may  be  highly  correlated,  and  the  error 
exponent  approaches  the  most  unfavorable  value  = 1. 
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Table  5 

ORDER  OF  LEADING  TERM  IN  ASYMPTOTIC 
EXPRESSIONS  OF  MODAL  OPERATIONAL  ERROR 


Integrator 

Type 

Comp. 

Path 

Small-stepsize, 
low-frequency 
range  Q « 1 

Large-step  size,  high-frequency 
range  Q » 1 

Strongly  A-stable 
operators 

Weakly  A-stable 
operators 

POE  AOE 

POE  AOE 

POE  AOE 

Pair  of  one- 
derivative 
operators  (7) 

(O') 

0(1)  0(0''') 

0(1)  0(0*^) 

0(w)  0(0^^") 

(0) 

0(1)  0(0'') 

(1.  2) 

o 

1 

0 

Two  -derivative 
operator  (Alb) 
(conventional 
implementation) 

(OM) 

o(q'^)  0(0'^"') 

0(1)  0{Q^) 

0(0)  0(0^^'') 

(2) 

0(1)  0(0^') 

NOTES: 

1.  POE:  Modal  amplification  factor  (43)  of  propagated  operational  error, 

2.  AOE:  Modal  amplification  factor  /Au  of  accumulated  operational  error 

where  AvT  is  an  statistical  average  of  local  modal  errors  j Au^[. 

3.  The  error  exponent  u varies  from  1 for  perfectly  correlated  sequential  local 
errors  to  1/2  for  uncorrelated  (random)  sequential  local  errors.  (The  value 
appropriate  to  the  low-frequency  range  is  not  necessarily  equal  to  the  value 
appropriate  to  the  high-frequency  range,  but  virtually  nothing  is  presently 
known  as  to  the  dependence  of  V on  Q.  ) 

4.  A strongly  (weakly)  A-stable  operator  is  characterized  by  the  largest  i oot 

of  the  polynomial  + . . . + having  a modulus  less  than  (equal  to) 

unity.  This  root  characterizes  the  filtering  properties  of  the  operator  as 
regards  suppression  of  high-frequency  noise  [4]  . 
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Section  7 

INHERENT  ERROR 


A Sample  Nonlinear  Equation 

As  previously  noted,  the  inherent  error  arises  from  uncertainties  in  the 
forcing  term.  As  a general  rule,  the  effect  of  inaccuracies  in  the  specifi- 
cation of  external  forces  can  be  ignored  in  practice,  unless  the  external 
force  field  results  from  the  concurrent  solution  of  a coupled  initial  value 
problem  (e.  g.  , fluid-structure  interaction,  coupled  thermoelasticity). 
Internal  force  errors,  on  the  other  hand,  may  be  important  in  nonlinear 
structural  dynamics,  as  the  following  sample  problem  illustrates.  Con- 
sider an  undamped  nonlinear  dynamic  system  governed  by  the  state 
equilibrium  equation 

M u + N(u)  = 

where  N(u)  is  a nonlinear  stiffness  term.  Two  widely  used  linearization 
techniques  for  treating  that  term  are 

N(u  ) = N(u  J K , (u  - u ,)  - AN 
i:'— n-1'  /-^n-1  '-m  — n-r  — n 

(68) 

— *“Ti  ‘ n *“  n ““n  n — i n ~ z.  n 


where  K is  the  Jacobian  of  N(u)  evaluated  at  u , K is  a Jacobian 
~n-l  — n-t  ~ 

evaluated  at  the  reference  state  u = 0_,  Q is  a nonlinear  "pseudo-force" 
term,  cp  is  a extrapolation  parameter  in  the  range  0 to  1,  and  is  the 
local  approximation  error.  The  Taylor  series  linearization  (68a)  is 
called  the  tangent-stiffness  method,  whereas  the  fixed-point  linearization 
(68b)  is  referred  to  as  the  pseudo-force  or  secant-stiffness  method  in  the 
engineering  literature.  For  the  sake  of  brevity,  we  use  the  latter  technique 
in  the  following  analysis.  Substitution  of  (68b)  into  (67)  evaluated  at  t^ 
yields 

Mvi  + Ku  = f + Af  (^>9) 

~ -n  -n 
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(70) 


where  the  effective  force  f and  inherent  error  source  ^f  are  given  by 

— n 

f = f +(l-hcp)Q  ,-cpQ  9 

Af  ==  Af  + AN 

— n — n 

In  (70b),  Aj^  denotes  the  applied  or  external  force  error,  as  in  (33b).  The 
practical  importance  of  the  term  AN^  should  not  be  overlooked.  In  problems 
involving  path -dependent  material  nonlinearities,  such  as  plasticity,  the 
extrapolation  error  in  (68b)  may  lead  to  a pseudo-force  error  of  the  order 

of  1-10%  of  f unless  a reliable  equilibrium-correcting  strategy  is  used. 

— n 

In  a highly  nonlinear  problem,  the  inherent  error  source  (70b)  propagates 
in  a complex  manner.  As  the  modal  spectrum  changes  with  time,  the 
principle  of  superposition  of  individual  sources  is  not  applicable.  The 
following  considerations  are  restricted  to  linear  systems,  but  should  also 
provide  insight  into  the  "semilocal"  force-error  propagation  in  mildly 
nonlinear  problems. 

Liinearized  Analysis 

Eq.  (69)  is  linearized  by  assuming  that  the  right-hand  side  terms  are 
independent  of  u.  The  propagated  inherent  error  q(t)  associated  with  a 
single  source  /\f^  at  t = 0 then  satisfies  the  linear  difference  equation 
(39b).  The  appropriate  initial  conditions  are 


tl 

^6  ^ 

% = 

(71) 

% = 

II 

1 

CTl 

3.-2  ^ - 

All  of  the  results  obtained  in  previous  sections  concerning  operational 
error  propagation  apply  to  the  inherent  error  function  ^(t)  if  and  p^ 
are  replaced  by  and  given  by  (71b,  c).  It  is  not  convenient,  however. 
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to  use  the  analog  of  (43)  as  modal  amplification  factor.  A more 

ITlaX  O 

significant  measure  is  the  ratio  of  modal  force  errors: 


I + a 


which  agrees  with  high-frequency  end  Q = uuh  » 1. 

Inherent  Error  of  the  Average  Acceleration  Method 

The  average  acceleration  (AA)  method,  which  is  a member  of  the  two- 
derivative  Newmark's  operator  family  [6]  , is  presently  the  most  commonly 
used  implicit  integration  scheme  in  linear  structural  dynamics.  It  is  also 
widely  used  in  nonlinear  dynamics  notwithstanding  its  poor  performance  as 
regards  high-frequency  noise  amplification,  experimentally  noted  by  many 
investigators  [9,14]  . The  error  growth  is  caused  by  the  coalescence  of 
the  two  characteristics  roots  of  the  integration  operator  as  Q = ujh  co 
if  either  computational  path  (0)  or  (1)  is  followed.  Near-root  coalescence 
brings  out  a "numerical  beating"  phenomenon,  which  is  briefly  analyzed 
below  and  substantiated  in  the  numerical  experiments. 

The  propagation  difference  equation  (39b)  associated  with  the  AA  method 
is  (cf.  Appendix  C ): 

% - + "In-a'  = » 

The  general  solution  of  (73)  is 

'In  = ^1  + ^2  ^2  (-74) 


where 


X,  ^ - e , 0 = tan  = 

1-1/4Q^ 


The  appropriate  initial  conditions  for  (73)  are  q = 0,  q = q , whence 

-1  o o 

(74)  becomes 

•■n  = % I - '■'2  ‘ ' ' <■ 
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which,  upon  introducing  (75),  can  be  reduced  to 

q = q Q^(l  + l/4  Q^)  sin  n 0 
^n  ^o 

The  beating  frequency  0,  given  by  (7  5),  approaches  4 Q ^ if  0 » 1 . 
amplification  factor  (72)  associated  with  (77)  is 

= Q (1  + 1/4  Q^)/(l  + n^)  1/4  fi  if  Q » 1 

The  same  beating  phenomenon  occurs  in  the  use  of  the  trapezoidal  rule 
when  implemented  in  the  (CO^)  form.  In  this  case,  the  appropriate  initial 
conditions  are  i ~ -q^  (since  h^  = l/2h),  which  causes  the  propagation 
function  (77)  to  double. 

Identical  conclusions  apply  to  the  propagation  of  operational  error. 
However,  this  source  is  not  usually  important  in  the  high-frequency  range. 


(77) 
The 

(78) 
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Section  8 

NUMERICAL  EXPERIMENTS 


In  formal  paper  presentation,  numerical  results,  if  any,  are  generally 
offered  as  verification  of  the  correctness  of  the  analytical  derivations. 

The  implication  is  that  the  analytical  framework  was  somewhat  put 
together  magically  by  the  investigators,  who  then  walk  to  the  computer  to 
reap  the  deserved  confirmation  of  their  brain  labors.  As  detailed  in 
Appendix  D,  however,  computationally -oriented  research  seldom  conforms 
to  such  a neatly  staged  schedule.  In  real  life,  numerical  experiments  and 
theory  development  proceed  concurrently,  with  the  former  often  suggesting 
avenues  of  progress  for  the  latter.  The  synergistic  interaction  of  experi- 
ment and  theory  has  in  fact  been  largely  responsible  for  the  success  of  the 
present  study.  With  these  sobering  thoughts  out  of  the  way,  we  proceed  to 
describe  summarily  the  procedures  followed  in  the  experimental  phase; 
this  is  followed  by  a detailed  discussion  of  the  numerical  results  presented 
in  Figures  3 through  27. 

Experimental  Setup 

The  numerical  results  presented  here  were  obtained  by  solving  the  propaga- 
tion difference  equations  (39),  specialized  to  an  undamped  linear  oscillator, 
at  a prescribed  set  of  sampling  frequencies  ranging  from  Q - a;h  - 10  to 
10  . The  study  parameters  considered  were:  (a)  the  integration  method, 

(b)  the  error  source  type  (operational/ inherent),  and  (c)  the  computational 
path  followed.  Observed  peak  responses  were  collected  into  amplification 
factor  arrays  for  subsequent  generation  of  amplification  spectrum  diagrams. 
In  addition,  options  were  provided  to  produce  time  history  plots  of  the 
propagatio*"  process  at  specific  sampling  frequencies. 

The  computer  program  that  carries  out  those  tasks  was  not  assembled  as 
a final  product,  but  evolved  in  stages  closely  integrated  with  progress  in 
the  theoretical  understanding  of  the  error  propagation  process.  The  initial 
version  was  used  to  investigate  the  error  propagation  characteristics  of 
one-step,  one-derivative  methods  under  the  basic  computational  paths  (0-2). 
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Multiple  step  one-derivative  operators  were  incorporated  once  the  basic 
properties  of  one-step  operators  were  reasonably  well  understood,  and, 
finally,  a set  of  two-derivative  operators  was  implemented.  The  distinc- 
tion between  paths  (0)  and  (O')  was  detected  in  the  final  stages  of  this  study 
(cf.  Appendix  D). 

The  present  version  contains  eighteen  one-  and  two -derivative  operators 
as  built-in  methods.  These  integrators  are  simply  accessed  by  symbolic 
key,  e.  g.  , METHOD  TR'  calls  for  the  trapezoidal  rule  coefficients.  (Some 
of  the  built-in  operators  are  actually  integrator  families,  such  as  Newmark's 
and  Wilson's,  in  which  case  the  identification  key  is  accompanied  h)  one  or 
more  parameter  values.  ) The  program  also  allows  the  direct  specification 
of  integrator  coefficients. 

i 

The  error  analysis  code  can  be  operated  in  either  batch  or  conver sational- 

interactive  mode  through  a free-field,  problem -oriented  language.  The  i 

interactive  capability  has  been  found  to  be  valuable  when  the  program  is 

used  as  an  experimental  tool  to  guide  and/or  confirm  theoretical  developments. 

Results  can  be  displayed  in  either  alphanumeric  or  graphs  mode.  Graphic 
display  is  generated  through  the  DISSPLA  plotting  par  ..age.  Graphic  output 
produced  in  interactive  work  can  be  directly  displayed  on  in-line  Tektronic 
terminals.  These  "quickie"  plots  are  useful  for  program  and  input  data 
debugging,  and  also  represent  a convenient  means  of  preparing  the  genera- 
tion of  high  quality  off-line  plots. 

Organization  of  Numerical  Results 

A comprehensive  set  of  numerical  results  on  error  propagation  characteris- 
tics of  specific  integration  methods  is  presented  in  PTgures  3 to  27.  Two 
set -inclusion  criteria  were  used:  (a)  the  result  pertains  to  a commonly  used 

integration  method,  or  (b)  the  result  illustrates  a particularly  important  or 
noteworthy  feature  previously  discussed  in  this  part.  All  of  the  plots  were 
entirely  computer -produced  on  the  Comp80  microfilm  recorder  through 
the  DISSPLA  plotting  package;  only  the  figure  captions  have  been  added. 

The  plot  series  is  organized  as  follows. 

Figures  3-11:  Amplification  spectra  of  propagated  operational  error  for 

one-derivative  methods 
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Figures 

12-17: 

Amplification  spectra  of  propagated  operational  error  for 
two -derivative  methods 

Figures 

18-19: 

Time  history  plots  illustrating  the  propagation  of  small- 
and  large-stepsize  operational  error  sources 

Figures 

20-22: 

Amplification  spectra  of  propagated  inherent  error  for 
one -derivative  methods 

Figures 

23-24: 

Amplification  spectra  of  propagated  inherent  error  for 
two-derivative  methods 

Figures 

25-27: 

Time  history  plots  of  propagation  of  inherent  error 
sources  illustrating  the  formation  of  "beats"  in  the  averag 
acceleration  method 

In  all  results  pertaining  to  one-derivative  methods,  the  same  operator  was 
used  in  the  discretization  of  u and  v. 


Operational  Error  Propagation  for  One-Derivative  Methods 


Figures  3 through  6 show  amplification  factors  corresponding  to  the  general 
one-derivative,  one-step  A-stable  integrator  family 


^ ^ 1 

— n-1 


h [f  u 
— n 


+ (1 


f)  i ,] 

—n-1 


1 ^ f ^1/2 


(79) 


for  the  values  f = 1,  0,  0.  60,  0.  55  and  0.  50.  The  limit  cases  f = 1 and  f = 

0,  5 correspond  to  the  backward  Euler  method  and  trapezoidal  rule,  respect- 
ively. These  four  plots  provide  a fairly  good  picture  of  the  transition  from 
a highly  numerically -damped  BD  operator  (backward  Euler)  to  a weakly  A- 
stable  operator  (trapezoidal  rule).  The  gradual  development  of  resonant 
amplification  in  the  large  stepsize  (high  frequency)  range  Q > 2 for  path 
(O')  is  especially  noteworthy. 

Figures  7 to  8 illustrate  two-step  methods.  Gear*s  two-step  operator, 
being  an  A-stable,  highly  damped  BD  scheme,  displays  amplification 
behavior  similar  to  that  of  backward  Euler.  Figure  8 shows  results  for 
the  two-step  member  (27b)  of  Park's  three-step  starting  family  (27). 


Figures  9 to  11  illustrate  three-step  methods.  The  unbounded  amplification 
of  Gear's  three-step  method  (Figure  9)  for  1/2  < Q < 2 is  typical  of  a 
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stiffly  stable  method,  for  which  the  instability  region  intercepts  a portion 
of  the  imaginary  frequency  axis.  Figure  10  shows  results  for  a method 
labeled  as  Formula  3 in  Jensen's  paper  [15]  ; the  original  stiffly-stable 
method  has  been  A-stabilized  by  addition  of  recommended  numerical 
damping  terms.  This  operator  is  interesting  in  that  it  is  the  only  one- 
derivative  method  for  which  path  (1)  has  been  found  to  be  inferior  to  path  (2) 
in  this  experimental  series.  Park's  three-step  method,  shown  in  Figure  11, 
displays  amplification  characteristics  similar  to  those  of  the  other  two  A- 
stable  BD  operators  (Figures  3 and  7)  despite  having  less  numerical  damping. 

In  all  cases,  an  excellent  agreement  with  the  predictions  of  the  asymptotic 
analysis  expansions  (51)  and  (62)  can  be  observed  in  the  small  stepsize 
region  Q < 0.  1.  In  fact,  the  agreement  is  surprisingly  good  even  in  the 
intermediate  stepsize  region  and  up  to  Q ^ 1 for  most  methods. 

Operational  Error  Propagation  for  Two-Derivative  Methods 

Error  amplification  plots  associated  with  two -derivative  implicit  operators 
show  at  most  two  curves  corresponding  to  paths  (0)  and  (2),  as  paths  (0), 

(O'),  and  (1)  coincide.  For  BD  operators  such  as  Houbolt's,  all  paths 
are  identical. 

Figure  12  shows  amplification  spectra  for  the  average  acceleration  method. 
This  operator  deserves  special  attention  because  it  is  the  most  widely  used 
integrator  in  linear  dynamics.  Moreover,  this  scheme  is  amenable  to  an 
exact  amplification  analysis,  typified  by  Eqs.  (73)  to  (77),  which  provides 
factors  valid  for  all  sampling  frequencies.  The  curve  for  path  (0)  has 
the  equation 

H = n (1  + 1/4  (80) 

P 

which  is  symmetric  about  the  Nyquist  frequency  Q = 2, 

The  central  difference  method  (Figure  13)  is  the  only  explicit  operator 
included  in  this  series  on  account  of  its  wide  application  in  shock  dynamics. 
(The  integration  of  the  propagation  equations  was  carried  out  through  a spe- 
cial branch  in  the  program,  and  the  label  'path  (0)'  has  in  fact  no  meaning.  ) 

It  is  interesting  to  note  the  smallness  of  the  amplification  factor  even  while 
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very  close  to  the  folding  frequency  Q = 2.  As  in  the  case  of  the  average 
acceleration  method,  an  exact  treatment  of  the  propagation  difference  equa- 
tions presents  no  difficulties  and  leads  to  the  expression 


H 

P 


Q 


-1 


(1  - 1/4 


2 -1/2 
n ) 


for  Q < 2 


(81) 


The  Houbolt  method  (Figure  14)  displays  the  excellent  high-frequency  noise 
suppression  behavior  that  can  be  expected  from  a BD  operator  with  high 
numerical  damping.  This  feature  has  contributed  to  the  popularity  of  this 
integrator  in  nonlinear  dynamics. 


The  linear  acceleration  method  (Figure  15)  is  a conditionally  stable  member 

of  tv/'o  operator  families:  Newmark's  (for  p = 1/6)  and  Wilson's  (for  6 = 1), 

1/2 

The  amplification  factor  becomes  infinity  at  the  folding  frequency  12 
3.464.  It  is  interesting  how  much  worse  path  (2)  is  for  this  particular 
method. 


Wilson's  one-parameter  operator  family  is  A-stable  for  0 > 1.  3660. 

Figures  16  and  17  illustrate  the  cases  9 = 1.  37  and  0 = 2,  0,  respectively. 

This  integrator  has  been  criticized  on  account  of  its  poor  error  propagation 

characteristics  when  implemented  as  a self-starting  method  [16]  . In  the 

usual  implementation,  the  solution  is  advanced  beyond  the  next  time  station 

(t  ) up  to  t + (0  - l)h  by  an  extrapolation  procedure;  the  method  then 

regresses  to  t through  interpolation.  High  error  amplifications  (e.  g.  , 

4 ^ 

R ^10  at  Q ^ 0.  1)  are  in  fact  common  to  all  of  the  so-called  "advanc ed" 

P 

LMS  methods.  In  our  study,  the  Wilson  operator  was  implemented  as  a 
regular  three-step  LMS  method  [11].  Figures  16  and  17  indicate  that  in 
such  a case  the  low-frequency  error  amplification  behavior  is  similar  to  that 
of  other  two -derivative  methods,  whereas  the  high-frequency  attenuation 
resembles  that  of  Houbolt's  method  (cf.  Figure  14),  The  numerical  results 
suggest  that  the  dependence  of  the  small-stepsize  amplification  factor  on 
the  parameter  0 is 


K 

P 


K 

P 


e + 0(1)  , 

6n"Ve^  + 0(1). 


for  paths  (0),  (O')  and  (1) 
for  paths  (2) 


(82) 


If  these  expressions  are  correct, 
cal  for  0 = 6^^^  = 1,  817. 


all  paths  would  be  asymptotically  icfenti- 


! 
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Time  History  Propagation  of  Operational  Error  Sources 

Figures  18  and  19  are  included  to  illustrate  the  temporal  propagation  of  a 
unit  operational  error  source  pertaining  to  small  stepsize  (Q  = 0,1)  and 
large  stepsize  (Q=  100)  normal  components,  respectively,  for  the  trape- 
zoidal rule.  The  propagation  histories  in  Figure  18  display  the  typical 
smooth  harmonic  behavior  of  low  frequency  components,  as  predicted  by 
Eqs,  (49)  and  (61),  The  transitional  propagation  period  (cf.  Figure  2)  is 
barely  perceptible  at  the  scale  of  the  plot,  since  it  lasts  only  one  step. 
Figure  19  illustrates  the  high  frequency  ’’numerical  beating"  behavior  asso- 
ciated with  path  (O');  histories  for  the  other  paths  have  been  suppressed  as 
they  would  be  undistinguishable  from  the  time  axis. 

4 

Inherent  Error 

The  amplification  of  inherent  error  sources,  as  measured  by  the  force 

* 

error  ratio  (72),  cannof  be  significant  in  the  small  stepsize  range  on  account 

2 

of  the  presence  of  the  factor  Q . With  the  exception  of  weakly  A-stable 
operators  implemented  under  certain  computational  paths,  all  methods  dis- 
play bounded  amplification  as  Q = co.  Because  of  this  similarity  of  behavior, 
only  five  representative  cases  are  presented. 

Figures  20  through  22  show  inherent  error  amplification  spectra  associated 
with  the  backward  Euler  scheme.  Park's  three-step  method,  and  the  trape- 
zoidal rule,  respectively.  These  plots  are  chosen  to  typify  the  behavior  of 
strongly  and  weakly  A-stable  one-derivative  methods.  The  only  unbounded 
growth  occurs  for  the  trapezoidal  rule  if  path  (O')  is  used,  in  which  case 
H ~ 1/2  n for  » 1. 

q 

Figu  es  23  and  24  depict  inherent  error  amplification  spectra  of  the  average 
acceleration  method  and  of  Wilson’s  operator  (9  = 1,  37),  respectively.  The 
former  exhibits  unbounded  growth  for  paths  (0),  (O')  and  (1),  as  predicted 
by  Eq.  (78). 

• Finally*,  Figures  25  through  27  are  included  to  illustrate  the  process  of 
formation  ®f  high  frequency  "beats"  for  the  average  acceleration  method. 

A unit  modal  force  error  is  introduced  at  t = 0,  At  a sampling  frequency 
fl  = uah  = 1 (~  6 samples/cycle),  the  propagation  histories  of  Figure  25  still 
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exhibit  rough  harmonic  behavior.  The  incipient  emergency  of  beats  in 
path  (0)  at  Q = 10  can  be  perceived  after  some  study  of  Figure  26.  The 
beats  are  perfectly  developed  at  Q = 100  (Figure  27). 
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! Section  9 

i 

I APPLICATIONS 

We  have  so  far  examined  the  effects  of  implementation  forms  on  the 
behavior  of  the  propagated  computational  error.  This  error  is  different 
from  the  algorithmic  error  as  commonly  assessed  in  terms  of  numerical 
^ damping  and  phase  distortion.  From  the  foregoing  analysis,  it  is  clear 

I that  the  computational  error  should  be  taken  into  account  in  the  design  of 

I reliable  transient  analysis  software.  As  this  aspect  has  rarely  been  dealt 

I with  in  present  structural  dynamics  programs,  we  sketch  here  how  the 

results  of  the  error  analysis  can  be  put  to  use  in  the  development  of  step- 
size  control  strategies.  We  consider  the  simplest  case  of  linear  struc- 
tural dynamics,  in  which  the  bulk  of  the  excitation  energy  is  associated 
with  low-frequency  components,  and  assume  that  the  (JO)  form  has  been 
implemented.  A rough  idea  of  the  peak  (relative)  computational  error  can 
be  obtained  from  (34)  and  Table  5: 

(83) 

where  u:  . is  an  estimate  of  the  lowest  frequency.  A stepsize  h is 

min 

initially  selected,  usually  on  the  basis  of  algorithmic  accuracy  require- 
ments, and  E is  formed.  An  order-of-magnitude  estimate  of  the  condition 
number  C(E)  can  be  economically  obtained  once  E is  factored  [IZ]  . (Note 
that  C(E)  is  a function  of  h,'  and  varies  from  C(M)  as  h — 0 to  C(K)  as 
h OD.  ) The  local  accuracy  level  e can  be  estimated  as  n^  10  , where 

d is  the  digital  significance  used.  Because  of  the  high  uncertainty  on  most 
of  the  quantities  in  (8  3)  - -particular ly  the  error  exponent- -the  frequency 
ID  . needs  to  be  estimated  only  within  an  order  of  magnitude.  If  (83) 
exceeds  a desired  accuracy  threshold,  a higher  arithmetic  precision 
level  is  likely  needed.  In  any  case,  (83)  provides  an  idea  of  the  minimum 
acceptable  stepsize.  If  there  is  a capability  for  adjustable  precision  work, 
a much  sharper  estimate  of  (83)  can  be  obtained  by  carrying  out  a short 
integration  in  two  precision  modes. 


— "max 


o (uo  . h)'*"  C(E)  € 

' min  g 
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Another  application  of  the  error  analysis  is  the  problem -adaptive  selec- 
tion of  optimal  implementation  forms.  For  example,  if  the  mechanical 
system  is  undamped,  the  use  of  a pair  of  BD  operators  in  conjunction 
with  the  (Cl)  form  would  be  desirable  in  terms  of  minimizing  the  computa- 
tional effort  (cf.  Table  3).  A prior  estimation  of  the  computational  error 
determines  whether  such  strategy  is  permissible;  if  not,  the  integration 
should  be  carried  out  with  the  somewhat  more  expensive  (JO)  form.  Other 
situations  can  be  similarly  treated. 
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Section  10 
CONCLUSIONS 


The  material  covered  in  Part  II  represents  an  initial  contribution  to  the 
study  of  computational  error  in  direct  time  integration  of  second-order 
ODE  systems.  Although  the  application  to  structural  dynamics  has 
been  emphasized,  the  conclusions  readily  extend  to  other  physical  pro- 
blems modeled  by  stiff-oscillatory  systems.  A summary  of  our  main 
findings  follows. 

Implementation  forms  are  found  to  have  a decisive  effect  on  the  error 
propagation  behavior  pertaining  to  low-frequency,  small-stepsize  normal 
components,  while  the  algorithmic  properties  of  the  integrator  play  a 
comparatively  minor  role.  Both  factors  are  equally  significant  in  high- 
frequency  error  propagation. 

We  have  shown  that  the  best  error  control  is  achieved  when  one-derivative 
LMS  methods  are  used  in  conjunction  with  the  path  (0)  implementation.  In 
this  case,  the  local  error  propagation  equation  approaches,  as  the  step- 
size  is  reduced,  an  integro-differential  form  that  is  insensitive  to  local 
velocity  perturbations.  This  favorable  behavior  cannot  be  duplicated  by 
any  conventional  implementation  of  a two-derivative  method. 

The  trapezoidal  rule  and  the  average  acceleration  method  are  often  con- 
sidered to  be  optimal  integration  operators  in  linear  dynamics  due  to 
their  algorithmic  properties:  no  artificial  damping  and  minimal  phase 
distortion  within  the  class  of  A-stable  methods.  Under  certain  computa- 
tional paths,  however,  these  methods  exhibit  a numerical  resonance 
phenomenon  that  amplifies  high-frequency  error  components.  This  effect 
seriously  undermines  the  viability  of  such  implementations  in  the  frequent 
cases  where  an  effective  filtering  of  high-frequency  noise  is  required. 

It  should  be  stressed  that  these  conclusions  are  of  practical  importance 
only  if  the  computational  error  propagation  can  significantly  affect  the 


I 


quality  of  the  computed  response.  This  condition  is  seldom  easy  to 
identify  before  the  problem  is  presented  to  the  computer.  It  is  hoped  that 
the  information  presented  here  will  be  useful  for  diagnosing  whether 
accuracy  problems  due  to  error  propagation  exist,  and,  if  so,  in  providing 
realistic  estimates  of  the  local  precision  level  required. 

In  our  opinion,  two  areas  deserve  further  investigation:  (a)  the  propagation 

of  inherent  errors  in  highly  nonlinear  problems,  in  which  a significant 
exchange  of  "error  energy"  among  spectral  components  can  take  place,  and 
(b)  the  establishment  of  statistical  properties  of  sequential  local  error 
occurrences  in  large-scale,  ill-conditioned  systems.  The  study  of  either 
subject  is  likely  to  require  extensive  numerical  experimentation. 
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Appendix  C 

NUMERIC  QUANTITIES  FOR  SPECIFIC  OPERATORS 
(1)  Integrator  Coefficients  for  Some  One -Derivative  Methods,  Eq,  (7a) 


Operator^ 
(Stability ' ) 

m 

Coefficients  , 

II 

o 

. 

3 

o/  & 

i = 0 

i = 1 

i = 2 

i = 3 

Backward  Euler  (A) 

1 

O' 

1 

-1 

e 

1 

0 

Gear  2-Step  (A) 

2 

O' 

1 

-4/3 

1/3 

p 

2/3 

0 

0 

Gear  3-Step  (S) 

3 

O' 

1 

-18/11 

9/11 

-2/11 

p 

6/11 

0 

0 

0 

Jensen  F3’^''‘'(A) 

3 

O' 

1 

-1.  92601 

1.  13841 

-0.  21240 

p 

0.  52503 

-0.  02916 

-0. 29085 

0.  081 36 

Park  2-Step ' ' ' (A) 

2 

O' 

1 

-1.2 

0.  2 

0.  6 

0 

0 

Park  3-Step  (A) 

3 

1 

-1.5 

0.  6 

-0. 1 

p 

0.  6 

0 

0 

0 

Trapezoidal  Rule  (WA) 

1 

O' 

1 

-1 

p 

0.5 

0.  5 

(A)  = A-stable,  (WA)  = Weakly  A -stable,  (S)  - Stiffly  stable. 

Formula  3 in  [15]  ; includes  recommended  numerical  damping  to  attain 
A -stability. 

Two-Step  member  of  starting  family  (27). 
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(Z)  Integrator  Coefficients  for  Some  Two -Derivative  Methods,  Eq.  (Alb)  i 


Operator 

(Stability''''*') 

m 

Coefficients  , 6^  j 

i = 0,  . . . 

m 

y / b 

i = 0 

i = 1 

i = 2 

Hi 

Average  Acceleration  (WA) 

2 

y 

1 

-2 

1 

(^Newmark  p = 1/4) 

6 

1/4 

1/2 

1/4 

, Central  Difference  (E) 

2 

'Y 

1 

-2 

1 

6 

0 

1 

0 

Houbolt  (A) 

3 

'Y 

1 

-2.  5 

2 

-0.  5 

6 

0.  5 

0 

0 

0 

Linear  Acceleration  (C) 

2 

y 

1 

-2 

1 

(^Newmark  p = 1/ 

1 

6 

1/6 

4/6 

1/6 

Wilson,  6 = 1,  37  (A) 

3 

y 

1 

-2. 27007 

1.  5401 5 

-0. 27007 

6 

0.  31282 

0,  36820 

0,  05507 

-0.  00616 

Wilson,  9 = Z (A) 

3 

*Y 

1 

-2.  5 

2 

-0.  5 

6 

2/3 

-5/12 

1/3 

-1/12 

* ^ the  notation  of  Appendix  A. 

(A)  = A-stable,  (WA)  = Weakly  A-stable,  (C)  = Conditionally  stable, 

(E)  = Explicit 

Also  coincides  with  Wilson,  9 = 1. 
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(3)  Coefficients  p and  a Occurring  in  Asymptotic  Expansions  (50)  and  (62), 
for  Specific  Operators  and  Computational  Paths 


Operator 

p / a'*' 

Path  1 

1 

(0) 

(O') 

(1) 

(2) 

Backward  Euler 

P 

0 

0 

1.  0 

1.0 

a 

1.0 

1.0 

Gear  2-Step 

D 

0 

0 

2.  25 

2.  25 

a 

1.  5 

1 . 5 

Gear  3 -Step 

p 

0 

0 

3.  36 

3.  36 

a 

1 .89 

1.89 

Jensen  F3 

P 

0 

0 

6.6  5 

3.  63 

a 

1. 905 

3,47 

Park  2-Step 

P 

0 

0 

2.  08 

2.  78 

a 

1.67 

1. 25 

Park  3-Step 

P 

0 

0 

2.  78 

2.  78 

a 

1 , 67 

1. 67 

Trapezoidal  Rule 

P 

0 

0 

2.  0 

4.  0 

0 

2.  0 

1 . 0 

Average  Acceleration 

P 

1.0 

4.  0 

Central  Difference 

P 

1.0 

Houbolt 

P 

2.  0 

2.  0 

Linear  Acceleration 

P 

1.0 

6.  0 

Wilson,  0 = 1.  37 

P 

1 , 37 

3.  20 

Wilson,  6 = 2 

P 

‘ ■ - — ■ 

2.  0 

1 . 50 

=1'  The  value  of  o is  given  only  if  p = 0,  i.  e.  , expansion  (62)  holds. 
Computational  path  distinction  does  not  apply  to  explicit  methods. 
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Appendix  D 

POST  SCRIPTUM:  THE  WAY  IT  WAS 

The  organization  of  the  main  body  of  this  paper  follows  the  traditional  style 
of  scientific  writing,  viz.  , the  presentation  of  a chiefly  analytical  body  of 
theory,  followed  by  its  experimental  substantiation.  Experienced  researchers 
know,  however,  that  the  processes  of  experimentation  and  theory  development 
usually  interweave  in  a complex  manner.  The  formal  presentation  of  com- 
putationally oriented  papers  is,  in  many  respects,  a holdover  from  the  pre- 
computer era.  Numerical  experiments  weie  then  performed  either  to 
demonstrate  the  validity  of  a proposed  theory  or  merely  to  tabulate  derived 
results.  Computational  mechanics  has  by  now  evolved  to  a stage  in  which 
exploratory  computations  cannot  only  offer  direct  insight  into  physical  pro- 
cesses (numerical  simulation),  but  may  lead  the  way  into  the  establishment 
of  mathematical  theories.  In  this  regard,  computational  Work  is  gradually 
complementing,  if  not  actually  replacing,  conventional  experimental  methods. 
Note  that  we  have  labeled  Section  8 of  Part  II  "Numerical  Experiments" 
rather  than  "Numerical  Results."  In  this  Appendix,  we  shall  attempt  to 
describe,  in  rough  chronological  order,  what  we  were  trying  to  find  out, 
what  we  did,  and  why,  and  when.  We  hope  that  the  readers  (and  the  writers 
will  become  such  after  this  is  written!)  will  identify  with  and  benefit  from 
this  candid  reconstruction  of  the  course  (and  curses)  of  a research  pro- 
ject. 

The  Setting 

As  stated  in  the  Introduction  to  Part  II,  our  main  objective  was  to  find  a 
computer  implementation  form,  within  the  class  of  one -derivative 
implicit  integration  procedures,  which  requires  minimal  computational 
effort  per  step  and  yet  achieves  satisfactory  error  propagation  control. 

Three  implementations  were  initially  considered:  the  (f  - Ku)  form,  the 
(Mu)  form,  and  the  conventional  implementation  of  the  Newmark  family; 
these  were  later  classified  as  (JO),  (Cl),  and  (C2)  forms,  respectively, 
in  the  nomenclature  of  Part  I. 


I 
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Another  departure  point  was  an  awareness  of  the  existence  of  techniques 
for  assessing  the  propagation  of  rounding  errors  in  the  conventional  numeri- 
cal solution  of  first-order  initial  value  problems.  Virtually  all  of  the 
existing  theory  has  been  developed  by  Henrici  [1,  2]  . A strenuous  journey 
through  Reference  [1]  revealed,  however,  that  the  class  of  computational 
errors  we  had  in  mind  was  not  dealt  with  at  all.  The  error  sources  we 
were  mostly  concerned  about  included  starting  errors,  uncertainties  in  the 
evaluation  of  mass  and  stiffness  matrices,  and  algebraic  solution  errors 
arising  from  inaccurate  matrix  factorizations  (in  linear  problems)  or  extra 
polation/iterative  schemes  (in  nonlinear  problems),  in  addition  to  the 
rounding  errors.  Moreover,  Henrici  did  not  stress  the  distinction  between 
propagated  and  accumulated  error,  neither  did  mention  the  frequency - 
decomposition  approach  on  which  most  of  our  work  was  to  be  based. 

With  these  "initial  supplies"  behind  us,  we  plunged  into  the  research 
labyrinth,  through  which  we  progressed  by  a combination  of  iterative 
experimentation,  theory  development,  and  some  degree  of  luck. 

The  Initial  Adventure 


As  a starting  point,  the  one-step  method  family 

u =u  +h[fu  +(l-f)u  ] 0<f^l 

n n-1  n n-1 


(Dl) 


was  used  to  treat  a single,  undamped  linear  oscillator 


u + u)  u = 0 
n n 


(D2) 


adopted  as  model  error  propagation  equation  (u  here  denotes  displacement 
error).  The  initial  conditions  were; 


u =1,  u,  = ...=u  =0 

O ' -1  -CD 

= l/(hf),  u_^  = . . . = u_^  = 0 


(D3) 


These  conditions  were  justified  as  follows.  An  assumed  unit  error  source 
in  the  displacement  occurs  at  a typical  time  station,  which  can  be  taken  to 
be  t = 0 for  conveniency.  To  study  the  propagation  of  this  error  source, 
we  simply  ignore  error  sources  at  other  stations  as  well  as  the  computed 


11-4  5 


r 


L 


solution  on  account  of  the  linearity  of  Eq*  (D2);  this  scenario  leads  to  (D3). 

Proceeding  through  the  first  time  step  of  the  (Cl)  form  (see  Table  1)  gives 

2-2 
a = -u)  u - -UD 

0 o 

h,'"  = i + (l-f)hu  = [1  - f(l-f)  uu^h^]/{h  f) 

1 o o 

h,^  = u + (1  * f)  h u = l/f 

1 o o 

= 1/f  i [1  - f(l-f)  h^]  (D4) 

u^  = g^/ (1  + f^  0)^  h^)  -*l  + l/f  asujh-*0 

~ (^1)  l/f  as  uoh  -»  0 

Another  cycle  produces  the  values 
u.,  ^ 1 -f  2/f  ) 

^ J as  ujh  “•  0 

^2  ~ ) 

The  preceding  calculations  suggested  that  the  displacement  error  source  J 

(D3a)  is  amplified  on  subsequent  stations  if  the  (Cl)  form  is  used.  The  ^ 

next  obvious  question  is:  how  does  the  peak  amplitude  depend  on  uoh?  j 

This  question  was  posed  to  the  computer.  In  this  instance,  numerical  | 

experiments  were  preferred  over  analytical  procedures  on  account  of  |] 

expediency  (quick  answers  are  very  important  in  the  initial  stages  of  a j 

research  project),  and  the  likelihood  of  crippling  mistakes  being  committed  j 

in  long  hand  calculations.  The  computer  output  indicated  that  peak  ampli- 
tudes were  proportional  to  (ouh)  in  the  small  stepsize  range,  a result 

obviously  in  harmony  with  the  i*  ^nd  suggested  by  (D4)  and  (D5),  j 


(D5)  j 

1 


The  (JO)  form  was  then  programmed  and  run  through  the  amplification  j 

analysis.  The  results  were  puzzling.  Peak  amplitudes  for  the  trapezoidal  | 

1 

rule  (f  = 0.  5)  and  backward  Euler  method  (f  = 1)  remained  very  close  to  two  i 

and  one,  respectively,  even  as  uoh  -*  0.  Back  to  pencil  and  paper.  One  ^ 

cycle  through  the  (JO)  computational  sequence  (Table  2)  gave  I 


l/f  I 

] as  o)  h 


0 


(D6) 
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and  the  same  asymptotic  values  recurred  cycle  after  cycle.  With  these 
preliminary  findings  in  our  possession*  it  was  time  to  begin  establishing 
a more  comprehensive  theory. 

The  Reasoning 


The  error  propagation  analysis  procedure  typified  by  Eqs.  (D4)  through 
(D6)  relies  on  the  step-by-step  integration  of  the  difference  equations,  after 
which  ojh  is  made  to  approach  zero  in  order  to  isolate  leading  expansion 
terms.  This  approach  is  straightforward  but  laborious,  and  provides 
information  only  over  as  many  time  steps  as  the  analyst's  endurance  can 
withstand.  Moreover,  specific  operators  have  to  be  introduced  ab  initio, 
which  is  hardly  conducive  to  the  establishment  of  general  conclusions.  On 
the  other  hand,  if  the  propagation  response  could  be  asymptotically  charac- 
terized by  the  closed  form  solution  of  a differential  equation,  not  only  would 
the  results  be  applicable  over  the  whole  time  domain,  but  the  tedious  case- 
by-case  approach  would  be  circumvented. 

We  were  therefore  motivated  to  seek  the  conversion  of  the  propagation 
difference  equations  to  a limit  differential  form  approached  as  ujh  -»  0, 
Propagation  histories  produced  by  the  initial  experiments  suggested  that 
u(t^)  approaches  a smooth  harmonic  behavior  (cf.  Figure  18).  Intuitively, 
one  would  expect  the  difference  equation  (15),  specialized  to  an  undamped 
oscillator,  to  approximate  the  model  O.  D,  E. 

u uo  u = 0 (D7) 

whose  solution  satisfying  the  initial  conditions  (D3)  is 

u(t)  = cos  ujt  + (ujfh)  ^ sin  uut  (D8) 

As  cjuh  -*  0,  the  sine  term  is  dominant  and  gives  a peak  of  order  (ujh)"\ 
which  is  the  correct  dependence.  This  result  is  not  particularly  startling, 
as  it  is  well  known  that  error  difference  equations  have  the  same  form  as 
the  governing  difference  equations.  Thus,  Eq.  (48)  was  easily  derived  as 
described  in  Section  5 of  Part  II. 

The  error  arresting  properties  of  the  (JO)  form  could  not  be  explained, 
however,  by  Eq.  (48).  The  numerical  results  indicated  that  the  solution 
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of  the  asymptotic  error  differential  equation  had  to  possess  the  structure 


u(t)  = u cos  uJt  + c sin  o)t 
o 


(D9) 


where  c must  be  0(1)  for  uoh  « 1.  Therefore,  the  limit  form  had  to  be 
insensitive  to  the  initial  velocity  jump  (D3b),  Fortunately,  we  were  aware 
of  the  properties  of  the  stochastic  differential  equation  used  to  model  noise- 
free  regulators.  Optimal  noise-free  regulators  possess  the  important 
property  of  smoothing  white  Gaussian  noises,  which  are  formal  derivatives 
of  Brownian  motion  processes.  In  stochastic  filtering  theory,  such  behavior 
is  modeled  by  the  Wiener-Hopf  integro-differential  equation.  As  the  "error 
velocity"  jump  (D3b)  can  be  viewed  as  the  formal  derivative  of  the  initial 
displacement  error  source,  it  was  conjectured  that  the  limit  error  propagation 


equation  of  the  (JO)  form  had  to  be 

t 

jm 

2 


U + UQ 


■f 


dt  = 0 


(DIO) 


And  so  it  was.  Setting  down  all  of  the  details  of  the  formal  derivation  of 
Eq.  (60)  required  many  hours  of  brainwork,  however. 

Getting  It  All  Together 

After  .'•he  important  milestone  of  deriving  the  asymptotic  error  propagation 
equations  (48)  and  (60)  was  reached,  effort  was  directed  toward  the  identifi- 
cation and  categorization  of  local  error  sources.  This  finally  resulted  in 
Eqs.  (32)  and  (33).  An  appropriate  term  was  needed  to  identify  displace- 
ment error  sources  such  as  (D3a),  which  become  important  in  the  small 
stepsize  range.  After  some  discussion  the  descriptive  name  "operational 
error"  was  agreed  upon. 

By  this  time,  the  error  analysis  program  had  grown  from  a small  deck  to 
one  box  of  cards,  and  a set  of  multiple-step  LMS  operators  had  been  imple- 
mented in  addition  to  (Dl),  A comprehensive  test  series  was  run  using 
several  integrators  and  the  three  basic  implementations.  An  examination  of 
the  output  revealed  that  the  peak  responses,  although  evincing  the  correct 
power  dependence  upon  ujh,  did  not  match  the  results  of  the  asymptotic 
analysis  as  far  as  the  proportionality  coefficient  was  concerned.  A study 
of  the  initial  stages  of  the  propagation  history  clarified  this  discrepancy. 

The  initial  conditions  (D3)  may  be  regarded  as  a unit  impulse  or  "error 
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shock"  applied  at  t = 0.  In  computational  shock  dynamics,  it  is  well  knov^.m 
that  if  the  spatial  difference  operator  used  to  treat  the  numerical  shock  dis- 
persion spans  m grid  points,  then  post-shock  effects  will  be  felt  over  0(m) 
grid  points.  Similarly,  in  the  case  of  temporal  shock  dispersion,  it  takes 
0(Zm)  steps  for  the  error  shock  to  disperse  when  m-step  one-derivative 
methods  are  used  (the  factor  of  two  appears  because  we  are  dealing  with 
second-order  O.  D.E.  systems).  After  0(2m)  steps  following  the  initial 
impulse,  the  asymptotic  error  differential  equation  is  closely  satisfied  and 
a smooth  harmonic  response  is  observed.  When  the  values  (42)  were  used 
as  initial  conditions  to  be  matched  by  the  asymptotic  propagation  solutions 
(49)  and  (61)^  an  excellent  agreement  between  theory  and  experiments  was 
reached.  The  physically  suggestive  term  "transitional  propagation  period" 
was  coined  to  identify  the  "temporal  boundary  layer"  shock-dispersion 
region  that  bridges  the  impulse  (D3)  with  the  free  propagation  oscillatory 
regime.  We  remark  that  this  procedure  is  a particular  case  of  the  "method 
of  matched  asymptotic  expansions"  [I'^l- 

The  Beaten  Path 

The  computer  program  was  farther  expanded  by  the  addition  of  two-derivative 
operators  (which  turned  out  to  be  a disappointment  from  the  standpoint  of 
operational  error  propagation  control),  the  implementation  of  inherent  error 
analysis,  and  the  incorporation  of  on-line  graphic  display  capabilities. 

Eventually  came  the  time  to  examine  the  propagation  of  inherent  (force) 
errors.  Analytical  considerations  had  indicated  that  the  propagation  of  this 
type  of  error  sources  can  be  significant  only  for  large  uoh,  i.  e.  , in  the 
high  irequency  range.  The  sensitivity  of  the  average  acceleration  method--  ; 

as  implemented  by  Chan,  Cox,  and  Benfieid  [l8]--to  the  propagation  of  ; ^ 

high  frequency  error  sources  had  been  noted  by  other  investigators.  The 
appearance  of  the  associated  numerical  beating  phenomenon  was  therefore  ^ ■ 

not  surprising,  and  was  readily  explained  by  the  elegant  analytical  treatment 
given  in  Section  7 of  Part  II.  The  same  phenomenon  was  observed,  however,  j 

in  the  trapezoidal  rule  implemented  under  the  programmed  (JO)  form;  this 
was  unexpected  and  led  to  a careful  re-examination  of  the  implementation  | 

procedure.  * 

( 

j 

j 
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The  Fixing 


While  the  authors  were  pondering  over  the  trapezoidal  rule  beating  in  path 
(O)--supposedly  the  best  error  control  implementation --a  quick  fix  was 
worked  out:  choose  B = D + K in  (23b)  instead  of  B = D.  This  trick  was 
successful  in  that  it  eliminated  the  beating,  but  turned  out  not  to  be  needed 
at  all!  It  was  discovered  that  the  (JO)  form  had  actually  been  programmed 
as  follows: 


u 

n 


g 


n 


u 

n 


n n + 1 


2 

UL)  u 


n-1 

- UL)  u 

n 

-1 

n-1 

h""  + 

n-1 

h V 
6 n-1 

(Dll) 


At  this  point,  v , appears  to  be  the  equivalent  of  u . for  the  model  equa- 
^ n-1  n-i 

tion  (D2),  as  the  mass  is  unity.  In  the  interest  of  saving  a few  words  of 
memory,  the  program  developer  fell  into  the  trap  of  using  the  same  array 
to  store  u and  v,  which  implies  the  (uncoded)  operation 


n-i 


n-1 


(D12) 


The  "past  velocity  correction"  obviously  violates  the  (JO)  sequence  listed 
in  Table  2.  The  investigators  of  implementation  procedures  had  fallen  prey 
to  an  implementation  detail! 

Once  the  four  state  quantities  u,  u,  v,  and  v were  carefully  separated,  the 
(JO)  form  performed  as  expected  over  the  whole  ujh  range.  We  decided 
to  label  the  calculation  sequence  involving  (D12)  as  path  (O'),  which  occurs 
naturally  in  the  conventional  (C)  formulation.  By  then  Figure  1 had  been 
back  to  the  drawing  board  for  the  N-th  time,  and  we  had  four  path  offsprings 
instead  of  three,  and  we  realized  that  the  quick  (JO)  fix  had  actually  opened 
the  doors  to  another  computational  path  family  (cf.  Remark  in  Section  4 of 
Part  I).  As  in  all  mushrooming  undertakings,  there  is  an  appropriate 
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time  to  stop  and  appraise  what  has  been  accomplished. 


1 


I 

1 

i 

The  Ending 

This  study  was  initiated  in  earnest  in  August  of  1976  and  continued,  inter- 
spersed with  other  activities,  for  the  next  five  months.  By  December, 

Part  I had  been  rewritten  once,  and  Part  II  twice  over.  Several  potentially 

interesting  topics,  such  as  the  statistical  accumulation  of  propagated  error 

sources  and  the  inherent  error  propagation  in  highly  nonlinear  problems, 

had  been  hardly  touched  upon.  Nonetheless,  1976  was  rapidly  drawing  to  a 

close;  it  was  time  to  prepare  for  the  1977  ASME  meetings,  the  Independent  ^ 

Research  and  contract  reports  would  be  soon  due,  and  the  error  analysis  | 

program  hovered  close  to  4000  cards.  So  we  decided  to  sit  down  and  finalize 

the  formal  paper! 


AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  ; BACKWARD  EULER 


CIRCULAR  SAMPLING  FREQUENCY  cjh 


Amplification  spectra  of  propagated  operational  error 
for  the  backward  Euler  method  [Paths  (O')  and  (Z) 
coincide  with  (0)  and  (1),  respectively] 


AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  : GOD  1-STEP  . F = .600 
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Fig,  4 Amplification  spectra  of  propagated  operational  error 
for  the  general  one-derivative,  one-step  method  (79) 
with  f = 0,  60 
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Fig.  5 Amplification  spectra  of  propagated  operational  error 
for  the  general  one -derivative,  one-step  method  (79) 
with  f = 0.  55 
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AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  ; TRAPEZOIDAL  RULE 
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Fig.  6 Amplification  spectra  of  propagatt-*i  . - 
for  the  trapezoidal  rule 


AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  : GEAR  2-STEP 
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Amplification  spectra  of  propagated  operational  error 
for  Gear's  two-step  method  [Paths  (0*)  and  (Z)  coincide 
with  (0)  and  (1),  respectively] 


AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  : GEAR  3-STEP 
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Fig.  9 Amplification  spectra  of  propagated  operational  error 
for  Gear's  three-step  method  [Paths  (O')  and  (2) 
coincide  with  (0)  and  (1),  respectively] 


AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  ; JENSEN  FORMULA-3 
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Fig.  10  Amplification  spectra  of  propagated  operational  error 
for  Jensen's  Formula-3  [15] 
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AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  : AVERAGE  ACCELERATION 
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Amplification  spectra  of  propagated  operational  error 
for  the  average  acceleration  method  [Paths  (0),  (O'), 
and  (1)  coincide] 


AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  ; HOUBOLT 
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Fig,  14  Amplification  spectrum  of  propagated  operational  error 
for  Houbolt's  method  (all  paths  coincide) 


AMPLIFICATION  OF  PROPAGATED  ERROR  30URCE 
METHOD  : LINEAR  ACCELERATION 
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Fig.  15  Amplification  spectra  of  propagated  operational  error 
for  the  linear  acceleration  method  [paths  (0),  (O')  and 
(1)  coincide] 
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AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  : WILSON  , 6 = 1 . 370 


CIRCULAR  SAMPLING  FREQUENCY  wh 


Fig.  16  Amplification  spectra  of  propagated  operational  error 
for  Wilson's  method  with  0 = 1.  37  [Paths  (0),  (O')  and 
(1)  coincide] 


AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  : WILSON  . 0 = 2 . 000 


PROPAGATION  OF  UNIT  ERROR  SOURCE 
METHOD  : TRAPEZOIDAL  RULE 
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PROPAGATION  OF  UNIT  ERROR  SOURCE 
METHOD  : TRAPEZOIDAL  RULE 

wh  = 1 . 00+02 


AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  : BACKWARD  EULER 


Fig.  20  Amplification  spectra  of  propagated  inherent  error 
for  the  backward  Euler  method 


AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  : TRAPEZOIDAL  RULE 
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Fig.  21  Amplification  spectra  of  propagated  inherent  error 
for  the  trapezoidal  rule 


AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  : PARK  3-STEP 


AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  ; AVERAGE  ACCELERATION 


f I • I I « f f 9 W 9 9 W W f 

10‘‘  10®  10‘  10*  10* 

CIRCULAR  SAMPLING  FREQUENCY  wh 


Fig*  23  Amplification  spectra  of  propagated  inherent  error 
for  the  average  acceleration  method 


AMPLIFICATION  OF  PROPAGATED  ERROR  SOURCE 
METHOD  : WILSON  . 0 = 1 . 370 
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Fig.  24  Amplification  spectra  of  propagated  inherent  error 
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PROPAGATION  OF  UNIT  ERROR  SOURCE 
METHOD  : AVERAGE  ACCELERATION 

cyh  = 1.0 


Propagation  history  of  an  intermediate  stepsize  unit 
inherent  error  source  for  the  average  acceleration 
method 
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PROPAGATION  OF  UNIT  ERROR  SOURCE 
METHOD  : AVERAGE  ACCELERATION 

«h  = 10  . 
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Propagation  history  of  a large  stepsize  unit  inherent  error 
source  for  the  average  acceleration  method,  showing 
incipient  appearance  of  numerical  beats  on  path  (0) 
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